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АВЗТВАСТ 


This report is written to serve as a guide to those persons 
who, having no previous experience with Monte Carlo methods, wish to 
apply these methods to their own problems. Particular emphasis is 
given to techniques which are useful in dealing with problems concerned 
with the diffusion of particles (and gamma rays) in material media of 
some complexity, both from a geometrical and а nuclear standpoint. 
Included as an appendix are brief summaries of a variety of problems 
of the above-mentioned type to which the methods described herein have 


been applied successfully. 
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FOREWORD 


The present report is a summary of the Monte Carlo method as it 
applies to problems involving the interplay of neutrons and photons with 
bulk matter in geometric systems of varying complexity: It is intended 
to serve as an introduction and practical guide for the fast-growing 
group of people who are concerned with auch systems. 

A brief resumé of some of the unclassified problems successfully 
treated by the methods here outlined is included &s an appendix and may 
serve to emphasize the practicality and flexibility of these sampling 
procedures. А similar resumé of some of the classified problems which 
have been handled is issued separately as LAMS-2121. 

The gener&l method was originally developed by Fermi, Ulam and 
von Neumann; many others have contributed special techniques and devices. 
No attempt has been made to give all sources in the text, and naturally 


no claims to originality are made for the procedures included. 
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CHAPTER Т 


BASIC PRINCIPLES 


1. General nature of the problem. All problems treated in the 
present manual involve estimation of what percentage of particles emanat- 
ing froma given source, after undergoing specified processes in a 
material medium of known geometry, can be expected to terminate in certain 
stipulated categories. 

If all relevant probabilities are known for the elementary events 
in the "life history" of such a particle, the Monte Carlo method is 
applicable, and indeed is usually the only method available. 

Moreover, its technique is pre-eminently realistic, consisting in 
actually following each of a large number of particles from the source 
throughout its life history to its "death" in some one of the terminal 
categories, using the elementary probabilities at each stage of its 
career in determining its fate. 

The present state of development of high-speed digital computers 
permits the use of samples of a size sufficiently large to ensure satis- 


factory accuracy in most practical problems. 


ail- 
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2. Outline of procedure. In any particular problem, а particle 
is completely characterized by a set of parameters which are sufficient 
to determine its (probability) behavior in all situations it may encounter 
during its history. These always include its position and direction 
coordinates, and in most cases its energy. Much more will be said about 
these and other particle parameters in the sequel (the appendix to this 
report and LAMS-2121). 

Тре Monte Carlo method of dealing with problems of the kind we have 
indicated breaks up naturally into а well-defined set of subroutines, 
which we shall briefly describe here, postponing their detailed treatment 
to subsequent chapters. They аге schematized in the following generalized 
flow diagram (Fig. 1). This is only intended as & general guide to the 
chapters of the text, &nd is subject to many revisions, depending on the 
special circumstances of the problem. 

(c). The proper assignment of all particle parameters to a source 
particle involves the spatial and angular distribution of the source, and 
its energy distribution, in case of а non-monoenergetic source. 

(в). A special routine may be provided to determine the position 
of first collision in cases of high transmission, where it is desirable 
to distinguish between this and subsequent collisions. 

(В). А general routine is designed to decide whether a particle, 
starting with known parameters from an arbitrary point of departure, in 
а particular zone of the system either suffers a collision within this 
Zone or reaches the boundary of this zone on its line of flight without 


incident. The essential physical concept involved in this decision is 
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Source parameters 
(Ch. II) 


Parameters at forced 1st 
collision (Ch. III) 


Collision in or escape 
from zone decision (Ch. IV) 


Collision. Parameters at 
point of collision 


Escape. Parameters at 
point of escape 


Ome 


Escape from 
system 


[m ]-o 


Nature of collision 
(Ch. V, VI) 


Direction after 
scattering (Ch. УП) 


(5)— 


Classification of 


escapes (Ch. VIII) 
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the mean free path. The subsequent procedure depends оп the nature of 
this decision, as indicated in Fig. 1. The particle parameters at the 
point of collision or escape are computed before proceeding. In case 
of escape, one proceeds to (8) or to (є) according as to whether 

the escape is to ап adjacent zone, or from the system. 

(Y). The collision routine naturally depends upon the physics of 
the medium and is subject to the widest variation. Its objective is to 
decide the exact nature of the collision and the immediate fate of the 
particle after collision. This includes the eventuality that the particle 
may terminate its career under the conditions of the problem; if this 
occurs, a tally is made in the appropriate terminal category counter, 
and one returns, as in all cases of termination, to an entry (a) which 
leads to the source routine (c) for introducing a fresh particle. In 
the event that the particle is to be followed further, the essential 
information required before proceeding to (8) isa knowledge of the 
laboratory angle of deflection from the line of flight and the energy 
and "weight" of the particle after deflection. 

(8). A purely geometric routine determines the direction param- 
eters of the deflected particle, and leads to (В). At this point, all 
particle parameters should be evaluated as they obtain at the point of 
collision, after deflection. 

(€). In the event that escape from the system occurs in the (В) 


routine, terminal classification is made and one returns to (e). 


= 
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This is the general scheme of the method аз it operates in all 
problems we propose to discuss. The chapters that follow are designed 
to elaborate on each of these subroutines in detail, the necessary 
physical concepts and their preparation for computation being developed 
as they are needed. 

Before this, however, must come a brief discussion of random 
numbers and the "fundamental principle of Monte Carlo," upon which all 


else rests. 


3. Production of random numbers. It is necessary to have upon 
call some source of random numbers equidistributed on the interval 
OSr«1. Ideally, one might spin a wheel of uniform scale, and indeed 
there exist Lists +2) of such Hünbéra generated in truly random fashion. 
There are computational algorithms adapted to digital computers which 
appear to serve our purpose just as well. Опе сап find descriptions of 


such methods in the literature‘ 1:2) together with discussions of the 


(oe, articles in U. S. Department of Commerce , National Bureau of 


Standards, Applied Mathematics Series #12, Monte Carlo Method, 
Washington, D. С., 1951. 


2 
( те Rand Corporation, A Million Random Digite with with 100,000 Normal 
Deviates, Free Press Publishers, Glencoe, Illinois, 71955. 


D. H. Lehmer, Mathematical Methods in Large-Scale Computing Units, 


Proceedings of of & Second Symposium on n Large-Scale Digital Calculating | 
Machinery, Harvard University Press, Cambridge, Massachusetts, 1951. 


Herbert А. Meyer, ed., symposium on Monte Carlo Methods, John Wiley & 
Sons, Inc., New York, 1956. | 
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"tests of randomness" which have been applied to them. We do not enter 
here upon such questions since the method adopted will probably be 
dictated by the type of machine used and other practical considerations. 
A description of the method we have employed will be found in Ch. IX, 

8 3a, b. We should like to call attention especially to the short paper 
by von Neumann, cited in footnote (1), which contains some words of 
wisdom for those who may be troubled by the "state of sin" accompanying 


the use of deterministic "random numbers." 


1. The fundamental principle of Monte Carlo. Suppose that a 
homogeneous medium consists of nuclei of three different types A, B, and 
C, and one knows that, in the event of a colitefon of a neutron in the 
medium, the probability of collision with type A is .2, with B is 
З, and with C is the remaining .5. It is intuitively clear that, 
if а large number N of random numbers are produced, approximately 

.2 N will fall on the interval 03г<.2 

.3N will fall on the interval .2& г <.5 

.5 № will fall on the interval .5$r«1 
and that this approximation will improve with increasing N; indeed, 
this is the salient feature we demand of any scheme of random number 
production. It is clear therefore how reference to a random number can 


be used to decide which of the three types of nuclei is hit in the event 


of а collision. 
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One of the decisive features of digital computers now existing 
is their ability to make decisions at high speed, with no limitation on 
the number of logical possibilities involved. The typical Monte Carlo 
flow diagram is largely occupied with intricate decision features. 

As a first example of how a flow diagrami 3) schematizes the pro- 
cedure in the above example we may note Fig. 2. 

More generally, if Ejs see, E, are n independent, mutually 
exclusive events with probabilities Ру, see р, respectively, and 


Py E p, = l, we will agree that 
ру + 2а ё ue a <P, + ... ue 


determines E,. This is the "fundamental principle" insofar as it applies 
to the discrete case of а finite number of eventualities. 

We may use it to duci ah а heuristic approach to the continuous 
case, in which it is desired to determine from а random number one of а 
continuum of values of & variable x. In this connection it is conven- 
ient to bear in mind the fact that the latter case may always be regarded 
аз approximable by a large finite number of distinct cases; indeed, in | 


computation with a fixed number of digits, we are always actually 


concerned with the discrete approximation. 


| 
om all flow diagrams of this report we adopt the convention: x 


means х 20, x^ means x<0. A box containing r always denotes 
reference to the special subroutine generating the next random number 
of the sequence. 


2j7s 
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Fig. 2 


Suppose that we arbitrarily assign a variable x оп the interval 


ИЛ 


О = х сп to the events Еу, У E with the agreement that 


і - 1 © х <i represents the event E, . Let us construct & probability 


density function p(x) by the definition 


р(х) = Р; do b 
Thus р(х) will be а step ТРЕ like that of Fig. 3. Now 


suppose that we define the probability distribution function 


x 


(x) -f "OET 
о 
whose graph is a broken line as indicated in Fig. 4. Note that P(O) = О, 
P(n) = 1. Since Р(1) = p, +... + Pjs ме may interpret Р(х) to mean 
the probability of the inequality Ё $x, for х= і, i=l, ..., п. 


Moreover, it is clear that the equation 


x 
г = P(x) -f P(E) at 
O 
determines x uniquely as a function of r, in such a way that if r 
is uniformly distributed on the interval O $r «1, x falls with 
frequency p " in the interval і - 1 & х <i, thereby determining the 
event E, under our agreement. 


| 
me figures are drawn for the simple A, B, C example. 
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р(х) 


м осо A @ 


Fig. 3 


Р(х) 
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Ме may state at once the fundamental principle as it applies to 
the continuous case: If р(х) dx is the probability of х lying 


between x апа x + dx, with agx «b, and 


b 
| p(t) а = 1 


Ф 


then 


r= (x)= | p(t) at 


determines x uniquely as a function of г; moreover, if г is 
uniformly distributed оп O$r«1, then x falls with frequency 
p(x) dx in the interval (x, x + dx). 

As a completely trivial first example, consider the case of 
neutrons that are to be located uniformly on an interval а $ x <b. 
We have p(x) dx = дх/(ъ-а), Г p(t) dé = 1, and г = Р(х) = 
f p(t) аё = (x-a)/(b-a), M x = a+r (b-a) determines x аз а 


а, 
function of the random number г. 


2. Application of the principle. It may be noted that the 


equation 


г = P(x) -| ple) а 
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can be expected to give rise to difficult implicit problems, since x 
must be determined from r. At worst, some successive approximation 
routine can be provided for the solution of the implicit equation 
г = P(x) when P(x) is obtainable in closed analytic form. Such time- 
consuming processes can be obviated in various ways, a few of which we 
indicate here. 

The simplest method, applicable in all cases, even. when P(x) 
is known only in experimental tabular form, consists in subdividing the 
(a,b) interval, storing accurate values of P(x, ) = Р, for the end 
points x, = а, < ху я x. - b of the subintervals, and using the 
discrete method for determining the subinterval (x 4» x, ) on which 
x falls, together with an interpolation for the actual value of x. 
If i is the first value of the index for which г - Pi is negative, 
r being the current random number, we may determine x from one of the 


formulas 


р-р Оч - ху.) 2) 


(2) 


(3) 
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Тре linear interpolation in (1) distributes x uniformly on the 
interval (x, > x, ) and is strictly valid only when p(x) is a step 
function. For sufficiently small subdivisions, (2) or (3) may give 
better results at the cost of an additional square root and are appro- 
priate when P(x) is concave up or concave down, respectively. 

One may contrast the latter formulas (2) and (3) with that 


resulting from a linear assumption for p(x) on (х,_,, х,): 


P, -pP (x-x y? 
1 1-1 1-1 


X. 2 


Pep so p, (x а x, i) + me 
i 1-1 


3. 


which is more complicated when solved for x, requires additional 
storage of the Pi » and uses & trapezoid rule for the Р.. 
Very useful also is a device employed by von Neumann, especially 
when р(х) is readily computable and storage space is at a premium. 
This consists of "throwing" points (¢,n) uniformly into the rectangle 
bounded by the lines х = а, х=Ъ, у = 0, у= 1 and rejecting the 


points lying &bove the curve 
X 
y = p (x) = р(х) /мах р(х) 


X being assigned the value & whenever (&,n) fall below the curve. In 
many trials, the ratio of the number of points retained with $ between 
x and x + dx to the number of points retained altogether will be 


approximately the ratio of the areas 
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* о x b 
р дау f р (%)4& = дах f р(%)4& = p(x)dx 


The method is illustrated by the flow diagram of Fig. 5. 

Obviously the device in this form is impractical if the area under the 
curve y= р (x) is a small fraction of the enclosing rectangle.  How- 
ever, modifications involving other “ООС areas can obviously be 
devised. 

One may also retain only points (€,n) &bove the curve y - a Gos 
assigning the value Ё to x and adjusting the "weight" of the particle 
by а factor р(х) f (1 - p (M/C - р G9). | 

Finally, а Р of the two methods may be used, а first 
random number determining the interval (x, 4, x, ) by reference to the 


Pi and the von Neumann device then used on р(х) on this interval. 


The method is then accurate, and the efficiency high. 


©) 
[Hem Hen HE 
© 


Fig. 5 
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CHAPTER ТТ 
THE SOURCE ROUTINE 


1. Introduction. It is the purpose of the present chapter to 
describe how a seeped is initiated by the machine, how "print-outs" 
are automatically effected, and how the particles are ee ee 
source. 

Suppose that we let N denote the number of particles already 
processed, so that М = O at the start of the problem. After the 
machine has processed any given number N of particles it will contain 
in various counters the numbers М, of these М whose careers have 
terminated in a set of disjoint, all-inclusive categories С.. Thus 
the ratios N, /N constitute the output of the problem апа serve as 
estimates of the probabilities of & source particle terminating in the 


various categories C It is ordinarily desirable to print the cumula- 


i’ 

tive totals М; periodically during the course of the problem, вау 
* 

every N particles, so that convergence and "reasonableness" can be 


observed. 


Thus the beginning of & flow diagram usually resembles that in 
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Fig. 6. It will be noted that N' denotes the number of particles 
processed since the last print-out, being reset to zero after each print, 
whereas М cumulates during the entire problem. The computation may be 
stopped at any time at which the stability of the N, /N or other 
statistical considerations indicate that sufficient accuracy has been 
attained. Тһе (а) entry, as has been mentioned, is that point to 
which the machine returns after the particle it was following has 


terminated its career in some one of the categories Ci. 


Print N, 
all Nj 


Fig. 6 
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Тре (о) exit leads to that part of the flow diagram devoted to 
assigning particle parameters to a source particle. 

There are two types of storage involved in all machine problems. 
"Permanent storage" places are reserved for constants like 0, 1, №, 
which do not change their values during the course of the problen, 
while "dynamic storage" refers to storage positions in the machine which 
contain values of parameters like N, N', м, which are subject to 
change as the problem progresses. 

It is desirable to keep a record of all storage of both kinds 
introduced into the flow diagram as it is constructed so that none be 
overlooked, and so that some estimate of the "size" of the problem can 
be gained as one proceeds. Occasionally it becomes clear that the 
memory of the machine is being exceeded, &nd various devices must be 
introduced for reducing the permanent storage (e.g., by multiple 


storage) or the length of the code itself. 


2. Particle parameters. From & consideration of the physical 
and geometric features of the problem, one fixes upon & set of particle 
parameters whose values at any time suffice to completely characterize 
the particle. We proceed to discuss these in detail. 

We shall limit our examples to two types of coordinate systems 
for space and direction; namely, we shall either employ space 
coorüinates x, у, Z, together with direction cosines u = cos a, 


у = cos В, w= cos у for direction of flight, where a, В, 7 are 


=272 
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the angles made by the line of flight with the x, у, 2 axes, respective- 
ly, or (in cases of spherical symmetry only) we shall use the radial 
distance К of the particle from the center О, together with the 

cosine м = cos y of the angle у which the directed line of flight 
makes with the positively directed radius vector. These coordinates 

are indicated in Fig. 7. We may note that OS a, В, у & m and 

on this range the cosines assume all values on the range -l < ц, v, 

м & + 1 once and only once. It is also helpful to remember that the 


= 


direction coordinates (и, у, м) may be regarded аз defining a point on 
the unit sphere us + ue + Ta = 1 in direction space U, V, W. 
Considerable advantages attend the use of parameters R, м in 
cases where spherical symmetry warrants it. However, we have found 
that in more complicated geometries, for example in cylinders, even 
when symmetry obtains, the use of coordinates indicated by the symmetry 
is hardly worth while. The main reason for this is that a particle 
which proceeds from some point of departure to a new place changes its 
directional coordinates in the latter case. The x, у, 2; u,v, w system 
has the great advantage that u, v, w remains unchanged under linear 
displacements. 
All problems require the use of spatial and directional coordinates. 
Other parameters are dictated by the nature of the problem. | 
Most problems are energy dependent; that is to say, the physical 


processes involved have elementary probabilities which аге functions of 


the particle energy. In such cases we carry ап energy parameter Е 
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Fig. 7 
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and usually ап energy group index g. Cross sections are usually too 
complicated as functions of energy to be accurately fitted by simple 
formulas, so that in practice the whole range of energy involved in 
the problem is subdivided into suitable groups g = l, ..., G by lower 
bounds Е. > E 


1 2 
tabulated for these intervals. 


AT > Bo? and &ll necessary functions of energy are 


Moreover, systems encountered in practice are frequently non- 
homogeneous, being composed of zones у= Дт 95 + of varying densities 
or of different materials. This usually involves storing physical 
quantities as functions of. $ as well as of g. Permanent storage 
then contains numbers "gy with two independent indices, which 
necessitates an additional particle parameter y to indicate the zone 
presently occupied by the particle. 

In problems involving high transmission, capture, fission, (п-2п) 
reactions or other such features, it is desirable, although not 
necessary, to introduce а particle parameter W, called its weight, 
which is initially unity &t the source. То illustrate its use, consider 
а problem in which, upon collision of а neutron with а nucleus, there 
is а probability р of capture. We have the alternatives of (a) not 
employing weights, using а random number г in case of collision to 
determine whether capture occurs or not, by reference to the capture 
probability р, and if г < p, losing the neutron to а capture cate- 
gory, returning to (a); or (b) using a weight parameter М, tallying а 
weight pW in the capture category deterministically, and continuing 


$ 
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with а neutron of weight (1-p)W, which now scatters on the proper 
nucleus. In the latter case we lose no trajectories to capture and 
get a better picture of the capture distribution itself. 

Certain problems are concerned with the time a particle takes to 
travel from the source to its death, which calls for a parameter т 
giving the "age" of the particle at all phases of its life. A parameter 
у is employed for the number of collisions suffered by а particle in 
some problems, e.g., in those dealing with thick target corrections. 

We include for easy reference a list of these most frequently | 
used particle parameters, апа will &dhere throughout to the notation 


indicated here. 


Particle Parameters 


Space coordinates X,y,2 or В 


Direction coordin&tes u,eV,W or м 


Energy E 
Energy group index g 
Zone number ў 
Weight W 
Age T 


Number of collisions у 
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3. Remarks оп units. It is perhaps desirable to mention briefly 
the matter of units at this point. We use the centimeter-gram-second 
systems of units, with the following qualifications. (а) Neutron епег- 
gies are expressed in electron volts (ev), Kev (103 ev), or Mev (106 ev). 
For the uninitiated, voltage has dimensions of energy per unit charge; 
one ev is by definition the energy acquired by an electron which has 
dropped through a potential difference of one (practical unit) volt. 
Since the latter is 108 /с of the electrostatic unit of voltage, (5) 
and the charge on the electron is е = 4.8025 x 1071? esu, it follows 
that one ev is eV = 4.8025 x 10-10 x 108/с = 1.60203 x 10-12 erg, 

(b) The energy of a photon is measured in (dimensionless) units of 
mc”, where m, is the rest-mass of the electron. This is explained 
more fully in Ch. VI. 

We include here the formula for computation of the time A T in 
seconds for а neutron of energy Е Mev to traverse a distance of d см. 
Letting k= 1.60203 x 1076 ergs per Mev, we have in the non-relativistic 


(6) 


range of neutron energies 


Ek = = mv 


(5) c= 2,99776 x 1010 cm sect, the velocity of light in vacuo. 
| | : 
(6) me relativistic kinetic energy is (m-m, )c , where m - n, / V 1- в? 


and В = v/c. As an exercise опе may determine the relative error 
(1.14) in computing v from the two formulas when E = 1h Mev, the 
highest neutron energy involved in the present report, which is thus 
confined to the non-relativistic range of neutron energies. We do 
give the relativistic treatment of the Compton effect in Ch. VI, 
however. | 
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where m,(gm) is the mass of the neutron and v its speed in cm ае 


Thus 


у = y2 E k/m 


Now one "gram atomic mass" of any physical particle contains 
Avogadro's number А = 6.0228 x 1023 particles. Taking 1.00893 Рог the 


atomic mass of the neutron, we have 


у =k! VE(Mev) сш sec ^1 


where k' = 13.83 x 109 numerically. Thus a 1 Mev neutron travels about 


14 meters per microsecond (и sec). Ме have then 


AT =k"d/ VE вес 


for the transit time AT, where К" = 7.231 x 10710, 


4, Space coordinates for source particles. We have indicated in 
§ 1 of the present chapter how the machine is led from the "start" of 
the problem to the point (о) at which it sets up a source particle, 
or having finished processing a particle, it is returned to (ø) via 
(м), while in $2 we have discussed the various kinds of parameters 


as they are carried throughout the career of a particle. We are now 
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ready to consider how initial values are assigned to these parameters 
as a particle issues from the source. We hope to give enough examples 
to indicate the nature of the procedure. 

(a). Uniform source on an annulus of radii К< Ry. Determina- 


1 
tion of coordinates x,y (cf. Fig. 8). The probability density function 


p(R) is 27 R/ т (Ro - R^), so that the fundamental principle sets 


R 
2 2 2 2 
г = Р(В) - | p(R) ав = (в - R )/(R, - Ro) 
R 
о 
Having thus located the radius К, we note that р(ф) dé = аф /2m 


во that the next random number determines ф by 


$a 
ait el 0+ уе 


We have therefore the routine indicated in Fig. O. 

(b). Uniform source in а spherical shell. The final result is 
R= 2 R2 + r(R> - R2) if spherical symmetry admits use of the single 
space coordinate В. If x,y,z must be specified, we shall have 
х = Ru, у = Rv, z = Rw, where (u,v,w) is а point uniformly distributed 
on the unit sphere üt + s + we = 1. How such points may be obtained 
will be discussed in the next section (5a) on direction parameters 


u, V, W, where the same problem arises. 


(c). Parallel-beam source incident on lateral surface of right 


circular cylinder of radius R,, height Н, or on a sphere of radius В 


1? ic 
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Fig. 8 
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If the beam is taken in the positive y-direction (u = 0, v = 1, w = O), 
and the cylinder in the position shown in Fig. 9, one may use the 
routine there indicated. 

For a similar beam incident on a sphere of radius Ri, we have 
X = в. yr, y=- y (RI -х)=- в, Ус - г), 2 = О, assuming symmetry 
about the у axis. 

(а). Isotropic point source at distance d from right circular 
cylinder. Here the determination of the x, y, z coordinates of the 
point of entry to the cylinder is correlated with the direction param- 


eters u, v, w. We therefore postpone this problem to the next section. 


5. Direction coordinates for source particles. We next consider 
the problem of assigning direction cosines u, v, м, or simply м in 
the. spherically symmetric case, to source particles. If the latter аге 
liberated in a material medium, the directions may be expected to be 
drawn from an isotropic distribution, whereas particles emanating from 
а surface are naturally limited to half of direction space and usually 
are distributed in some non-isotropic distribution about the normal to 
the surface. 

(а). Isotropic source; u,v,w.direction cosines. The problem 
involved is tantamount to that of choosing a point (u,v,w) uniformly 
distributed on the unit sphere „2 + y + ut = 1. The element of area 


for this sphere in spherical coordinates У, @ is sin» d» dé = 


- dw дф , where ф is the longitude, and we have used у instead of 
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Fig. 9 
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the traditional 0 for the remaining angle since у has already been 
introduced for this angle by way of w= cos y (cf. Fig. T). Те 
probability density function p(w) is therefore given by p(w) dw = 


- 27 sin» dy /un = 5 dw. We may therefore first determine w by 


and @ subsequently by 


$ $ 
r= | мө) авн f E en 
-N 


-7 


Reference to Fig. lO then completes the argument for the isotropic 
Source; we need only remember that P = s + = y - ut). 
(b). The cosine distribution. i This refers to & point source 
emanating from & surface. If we agree that the outer normal to the 
surface at the point has the direction u = 0, у = О, w= 1, then the 
cosine distribution has by definition the probability density function 


p(w) » 2w, with w z O. Thus 


(T) 


For & discussion of the way in which the cosine distribution governs 
particles emanating from a surface one may refer to Е. М. Sears, An 
Introduction to Thermodynamics, the Kinetic Theory of Gases, and 
Statistical Mechanics, Addison-Wesley Publishing Co., Inc., Cambridge, 
Massachusetts, 1955. Cf. the chapter on kinetic theory. 
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Fig. 10 
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апа м = yr in this case. The flow РМЕТ. Fig. 10 may be used 
with the replacement of w= 2г - 1 by w= уг. 

(с). Isotropic and cosine sources in spherically symmetric 
systems. In case spherical symmetry indicates the use of coordinates 
R, w defined in Fig. 7, we may have to assign the direction parameter 
w for point sources of the above kinds. Such sources will be located 
on spherical surfaces, and in case of a cosine source, the latter will 
be relative to the normal to such a surface, i.e., to the radius vector. 
It is therefore clear that the formula w= er - 1 is valid for the 
isotropic case, while w = Vr applies to outwardly directed, and 
м = -yr to inwardly directed cosine sources. 


(а). Isotropic point source at distance d. from cylinder of 


radius R,, height H (cf. Fig. 11). This problem presents some 


2? 


features of interest. Clearly, if we proceed naively to assign direction 


parameters u, v, w to particles as they leave the source 8, using the 
method of (a) above, most particles will fail to hit the cylinder, the 

size of which may be greatly exaggerated in the figure. The physical 
problem is naturally concerned with questions relative to the number of 
incident particles, not with the problem (an interesting one, incidentally) 
of what solid angle is subtended at the source by the cylinder. Ном the 
latter question may be attacked by Monte Carlo we leave as an exercise. 

In this connection it is worth mentioning that many purely geometrical 
problems which present forbidding analytical difficulties are apt to 


arise in photon problems, and may be (many have been) successfully 
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Fig. 11 
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treated by the simple sampling methods which we are discussing. 

To return to our point source, we observe first that the only 
particles which can hit the cylinder are limited to the wedge defined 
by the two tangent planes to the cylinder through the source S. This 
means that directions are limited, іп и, v, м space at the source to 
those with longitude 4$ оп the range  - Ф. < $ © $5 where Фо = 
arc sin R,/d, , and ф is equidistributed on this range. Moreover, for 
all those particles issuing from the source between ф and ф + dọ, 
the direction cosine м = cos y is limited by the end points of the 
element of the cylinder determined by ø, and on this range w has а 
constant probability density function, since the element of surf&ce area 
on the unit sphere is -dw dø. 

Now from the figure, the half-chord c, - Е ё (a, sin e* |, 
so that the distance d 


Thus for this ф, the range 


2 2 
в (н/2)/ V |(н/2) +a]. The 


routine for setting up х, у, Z, ч, у, м at the point of entry there- 


о = 4, cos @ - Coe 


of w is - W м 5 м; where w 


$ 


fore appears as in Fig. 12, where $5 » с» 42 › М> аге given by the 


above formulas. The last box follows from the fact that 


2/4. = tan(m/2 - у) = cot» = cosy/siny = м/р. 


(е). General distribution in half of direction-space. Occasion- 
ally it is necessary to consider a point source with w z O having 
some experimentally determined distribution in м, which may be given 


in the form of a table. We may then tabulate P, = 0< Pi iua РІ 2l. 


= 1 > м ... > W. = О, where Р, is the probability of a cosine 


1? I 
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Fig. 12 


w Ж w,, and proceed as in Fig. 13a (cf. also formulas (2) and (3), 
Ch. I, 85). When various distributions are to be tried, it is preferable 
to run а number of different problems, each for a specific м; the re- 
sults may then be weighted to give terminal percentages for any desired 
source. 

(f). А prejudiced source. It is sometimes desirable to sample 
certain ranges of source directions more thoroughly than others. To 
illustrate, suppose that we have an isotropic source, with w uniformly 


distributed on -1& м < l, but that the position of a counter makes 
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more important those particles which leave the source оп the range 

-L& wg Ww < 0. We may then give source particles equal likelihood 
of starting on this range or its complement, provided we assign weights 
(7+1) апа (1 - F), respectively. In this way about half of all 
source particles, of relatively smaller weight, originate in the 
important cone, the total weight processed for М particles having 
expectation > N(W +1) + 3 N(1- T) = М. The procedure is schematized 


in Fig. 13b. 


6. Energy of source particles. In the case of energy dependent 
problems, one usually decides upon & set of energy ranges with lower 
bounds Е > E, > ass > Ес, with storage of physical quantities for 
each of these ranges, Ес being the lowest energy permitted to particles 
in the problem. Particles which by chance reach lower energies are 
relegated to a terminal category reserved for such losses. An index 
g = l, ..., G designates the group number. If the source is mono- 
energetic, one simply sets Eo —» E, gy 8g, that is, E and g аге 
assigned thevalues of the initial energy Ej and the index gg of the 
group in which this energy falls; g, may or may not be unity; one 
may wish to study the behavior of а series of sources of various initial 
energies. 

If the source particles аге not monoenergetic but аге chosen 
from some given energy distribution, one tabulates P, =O< ... < P 


G 


= 1 together with E. >... > Bo and proceeds exactly as indicated 
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in §5e above, reading g for i, E for w,, and E for ж. 


Т. Other source parameters. The other parameters mentioned 
in Ch. II, 82 ; if called for, have the following values at the source: 
(a) >. is set equal to the number of the zone into which source 
particles enter (which may depend on the values assigned to spatial and 
directional coordinates, and thus involve a decision routine); (b) the 
weight М is usually unity at the source; (c) age т = О initially; 


and (а) » = О for the number of collisions undergone. 


8. Source for а -type calculations. Іп an important class of 


problems, including those which are designed to determine the rate of 
growth of a neutron population in a fissionable material, the source 
consists of a specified number of neutrons Np» 1 in group g; 

zone >, and direction range i. We will illustrate for the case of & 
homogeneous sphere of radius Ry» where these ranges are determined by 


bounds 


Bo > E > ооо > Es 
О = Во < К. < eee < К, 
1 = Wo > Wy >... > мр = е1 
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The object of such a problem is the determination of the distri- 
bution N E yi in which this initial distribution results &t time 
39» 


AT, &nd then to use the output м as the input source distri- 


yi 
bution for the next cycle. А SES PN of the (a) and 
(c) routines is indicated in Fig. 14. The determination of E, В, 
and м оп their ranges may be achieved by appropriate interpolation, 
or, if preferred, by deterministic &ssignment of suitable mean values. 
Ihe initial v&lues of Nguoi are arbitrary in а -calculations, but 
are guessed as well &s possible to hasten convergence to the limit 
distribution. The E» ^p i refer to the category being processed, and 
must be distinguished from the neutron parameters 853; i, which are 
subject to change while the category from which the neutrons arose is 


being processed. The feedback of N' for М is subject to suitable 


renormalization to preserve а reasonable source size. 
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(Cf. Ch. IV, $9) 


Fig. 14 
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СНАРТЕВ ТТТ 


THE MEAN FREE PATH AND TRANSMISSION 


1. The cross section concept. The cross section o ofa 


stationary "target" particle for particles of energy E, relative to 

a given single process may be thought of as the area presented by the 
target, assumed stationary in the laboratory system, to a beam of 
(point) particles of this energy, relative to the laboratory system. 
If we regard a thin slab of material of area 9€ , thickness а А ; numer- 
ical density N (target particles per cm? ), traversed by а parallel 
beam of particles (of energy E) normal to а, the total area presented 
to the beam by target particles is oNa (ад ), assuming ай so small 
that no "shadowing" ЕЕ, The fraction of beam-particles undergoing 
the process in this volume should be ома (а ү ) / a. We, therefore, 


have the attenuation law 


dn = -n No(ah) 


for the number n of particles in the beam, and 


n= n, exp ( -NoÀ) 
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represents the number of particles remaining in the beam after traversing 
a distance Å in such a uniform infinite medium, no being the number of 
particles in the beam at Д = 0. Note that we are assuming no other 
competing processes exist. 


It is, therefore, supposed that 
Que | exp Creo | но (ak) 


is the probability for а first collision between баа А+ ak , and 


P(À) = nn (-No 4) No (af) = 1 - exp (-но 4) 


О 


is the corresponding probability distribution function for а first 


collision at distance 2 }. 


2. The mean free path. The average distance A to first collision 


is defined as the first moment of the function p(À), i.e., 


CO CO 
A= J "E269 a. f exp (-NoÀ) Nof(af) = 1/Мо 


and is called the mean free path for the process at this energy. 
It follows that the Monte Carlo determination of distance } from 


an arbitrary point of departure to first collision, assuming the medium 
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homogeneous and infinite must be 


г = Р(ў) = 1 - exp (-4,/А) 


or 


Д - -à fn (1 - r) 


Since 1 - г is equidistributed оп O < r <1 if г is, we may use 


simply 


= -Afnr 


In practice, we are almost always concerned with a complex of 
different processes, each presenting its own cross section to the beam. 
specifically, we may have to deal with a medium containing different 
types A, B, С, ..., of nuclei; moreover, each type A nucleus may have 
а variety of different types of cross section, for example, an elastic 
scattering cross section оА(е1.), an inelastic cross section o^ (1n.), 
а cross section c^ ( fiss.) for fission, c^ сар.) for capture, and so 
on. 

The sum of the cross sections c^ (e1.) + e^ (in.) +... Of all 
types for а particular nucleus А is called its total cross section 
e (*ot.). If the medium contains nuclei of types A, B, C, ... in 
numerical densities N,, Мы, Nas ..., respectively, the "total cross 


| B 
section" for the medium is defined to be 
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r= N, (40%. + М; eP(tot.) +... 


Reference to the preceding discussion makes it clear that the 
the simple No of that argument should be replaced by > in the general 


case, and the mean free path for the medium is, therefore, 


x 3/3 


The Monte Carlo method is always concerned with the distance ) 
from point of departure to collision, and only in case of collision, 
turns to a consideration of the nature of target hit, and the type of 
process involved. Thus it is always the mean free path A= 1/ $ that 
is used and never the free path for any of the individual processes. 

It must be remembered that cross sections are, in general, ае- 
pendent upon the energy of the particle (picturesquely, the size of the 
target depends on the speed of the arrow); thus we write v ^(el.), ea 


c e (tot.), t and A_ as functions of energy E, and, in practice, the 


E 
free path is usually tabulated as Ag? g = 1, ..., G, where g is the 
energy group index. 

Moreover, in systems consisting of zones of differing composition, 


Y 


we will have an additional zone index on the free path, thus А в" 
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3. Ап example. Consider the problem of determining A for neu- 
trons of energy E = 3 Mev in a medium of CH, of density ё = ‚92 gm em P. 
At this energy, C has ап elastic scattering cross section ce" (e1.) = 1.14 
barns (1 barn = 1072% ene) while H has a similar cross section оН(е1.) 
= 2.23 barns. No other processes are involved. 

The atomic weight of С is 12 and of H is 1, so that the molecular 
weight of CH, is 12 + 2(1) = 14. Та one gram molecular weight G of any 
compound are А = .6 x 102% molecules. The mass of one CH, molecule is, 
therefore, G/A gm. Since 1 en? of CH, has mass 6 gm, the number of 
24 


molecules of CH, in 1 cm? is Н = 8/(G/A) = вА/С = .0394 x 10 сш”. 


Hence, the numerical densities о? С апа Н are № = N, Мн = 2N and 
Е = №, e" (tot.) + Ng c" (tot.) = N(o°(e1.) + 2 о М(е1.)) = .221 cm ts 
Thus finally A= 1/2 = 4.52 cm. 

In problems involving many zones of the same ETT at different 
densities, it may be necessary to store only the basic nuclear constants, 
together with zone densities, and provide for the machine to compute its 


own A Y when needed, by reference to the stored quantities. For instance, 


in the preceding example, taking energy dependence into account, we should 


have 
a2 = 1/6 y^ 
where K, = с (o: * 207) 
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and ё are stored quantities. Such a procedure is time-consuming, 
especially when it necessitates computation of the probabilities for 
type of collision upon each collision. 

The concept of free path as we have formulated it applies to 
photons in their interaction with electrons and nuclei as well as to 
neutrons interacting with nuclei. A discussion of photons will be 
found in Chapter VI. Throughout the report, except in Chapters V and 
VI, we speak in a general way of "particles" which may be photons or 


neutrons. 


4. "Small" systems and transmission. Consider a homogeneous 
medium having mean free path A for particles of a monoenergetic 
source. If the distance from source S to the boundary of the system 
along a given direction is L, then, of a beam of N source particles 
leaving the source in this direction, N exp (- L/ A) will escape un- 
deterred. It is clear that if the dimensions of such a system are 
small compared to the free path A, most source particles will escape. 
This is especially undesirable if the assignment of source parameters 
is complicated, and,in any case,requires needlessly large sources to 
produce an effective sample. 

Now there is nothing to prevent us from regarding а single 
mathematical particle leaving the source in а given direction as 
representing a large set of W actual particles. Since the output 


of all problems is & set of ratios N,/N, where N is the total number 


«5l. 


Google 


of source particles, the initial value assigned to М may be taken as 
unity. 

We may then argue that, for a particle of weight W (=1) leaving 
the source in the situation described above, a partial weight W exp (-L/A) 
is transmitted, the latter weight being tallied in a category T reserved 
for total transmission (without collision). Then we determine a position 
of first collision on the interval O < A <L within the medium for the 
remaining particle of weight м (1 - ехр (-L/a)) according to the formula 
г = P(f)/P(L), where Р(() = 1 - exp (-(/л), as derived in a preceding | 
section. Solving for у, one obtains Î = -a Ап (: - г[1 - exp (-L/a) ] +. 

If the medium is non-homogeneous, but may be regarded as consisting 
of a number of zones, each homogeneous in itself, the portion of the line 


of flight lying within the medium may be decomposed into successive in- 


tervals of lengths L 


q^ La? where зу is the segment lying in zone Y 


with free path А. The transmission is then clearly t = exp (-1, Л 1) 


v 


* exp (-L,/A,) ... EXD (-L /^ ,) = exp NS +... + "2 The 
probability of а first collision at а distance «1 S. L from the point of 


origin is, therefore, l - exp (-»), where Р is defined by 


1 ¥-1 uo 9 
and 
L L - +... + L ) 
Р = E TREE | q- + V Du eec ad 
A А-1 му 
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Thus р is the number of free paths represented by Q. 
therefore may record the weight Wt as transmitted, and force a first 
collision of the weight W(1 - t) on the line of flight at a distance 


А < L by means of the formula 


= (1- e^)/( - t) 


Thus P = - În [1 - r(1 - t)] determines ^y by means of the inequalities 
L L L 
pee uus. m 
1 7-1 N ^7 


ала. Q by the equation 


5. The "forced first collision" routine. We illustrate the use 
of the device in two typical problems in this and the following section. 
Consider a parallel-beam monoenergetic source incident on the lateral 
surface of a cylindrical shell of radii ROC i, and height Н, composed 
of homogeneous material having free path A at this energy. We suppose 
that the source routine has already assigned to a source particle its 
parameters at entry, вау, ч = 1, v = 0, м = O, and x, y, 2 (in the manner 


indicated in Chapter II, 84c), together with E = E , g = g,» v = О, and 


о? 
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М = 1. The exit from (с) should then lead to the (8,) routine for 
forced first collision as indicated in Fig. 16, based on the geometric 


properties of Fig. 15. 


Fig. 15 
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The procedure naturally depends on whether or not the line of 
flight crosses the hole. This is the reason for the |y | - Ro decision. 
The distance L is the total path length of the line of flight lying 
within the medium. The transmission T6 is based on the free path for 
initial energy group Bo: The distance traversed within the medium from 
point of entry to point of collision is denoted by Q. Note that the 
exit leads to the collision routine (у) with all parameters as they 


exist at the point of collision, momentarily before impact. 


6. Remark on the device in spherical problems. Although we 


shall not include an example of the forced first collision device in 
a system treated with spherical coordinates, we should mention that, 
if it is used to determine a weight W and a radial distance а the 
point of first collision, (8) one should exit to an entry such as (у!) 
discussed in Chapter IV, 83, which sets up the direction w and 

(У) itself. 

As an example of the forced first collision method in a non- 
homogeneous medium, consider a parallel beam of particles directed 
vertically upward (u = 0, у = 0, w= 1) and incident on a sphere 

2 2 2 2 


x ty +2 = Ro subdivided into spherical homogeneous shells by 


the radii В, > Ri Tb Ry 


| pr R,w are the source parameters, and À the distance from source 


to forced first collision, R2 = RÊ + 42 + 2RÍw in the solid homo- 
geneous case. A 


= O, zone pi having total mean free path 
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Fig. 16 


Nay For the incident energy. А set of storage places Pa s f, P, ате 
reserved in the machine for the numbers of free paths represented by 
segments (z',z") of the line of flight in the zones through which it 


passes. The method is illustrated in Fig. 16a. It is assumed that 


х = Ro vr, у = О have already been set for the point of entry in the 


source routine, symmetry obtaining about the z-axis. 
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Ap = Я 


Fig. 16a 
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Ap = (z" — z')/A 
A 


7. The transmission in subsequent history. The device of 


forcing collisions of the non-transmitted weight may, of course, be 
applied to collisions after the first. If applied consistently to 
all collisions, then no trajectory ever terminates in escape, and one 
would ordinarily rely on a weight cutoff for termination. Usually the 
geometric complexity of paths after first collision renders use of the 
device for further collisions impractical. We do not consider forced 


collisions other than the first in the present manual. 


8. Prejudiced first collision in "large" systems. Іп problems 
concerned with "large" systems, such as those arising in shielding, 
the transmission must be very small and yet one may have to obtain 
energy-angle distribution of escape and space distribution of various 
types of collision (e.g., inelastic collision апа radiative capture of 
neutrons in determining 7 -sources) throughout the system. Moreover, 
the existence of energy cutoffs makes very unlikely the arrival of 
particles in the farther reaches of the system after many collisions. 
In such cases, one may overcome the dwindling of first collisions due 
to the exponential by prejudicing the distribution of first collisions 
and weighting accordingly. 

We illustrate with а simplified example. Consider а thick plane 
slab of two layers (free paths Ais An at incident energy) bounded by 
the planes z = О, 2 = Lz = Lo» with a source directed vertically 


upward and incident on the surface z = О. Опе may then determine the 
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position 2 of first collision, together with the weight at this point 


by the scheme of Fig. 16b. Note that the formula z = rL, distributes 


first collisions uniformly throughout the slab, while the weighting 
gives the correct expectation for first collision between 2 and 2 + dz, 
namely, 


dz 2 =P -P dz 
e =е — 


2 ^n Am 


If N source particles are processed, the expected total weight 


assigned to these at first collision is then N(1 - T)» where T 
L La- L 
l 2 L NE 
не is the transmission. 
1 i 


Fig. 16b 
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СНАРТЕВ ТУ 


THE COLLISION OR ESCAPE ROUTINE 


1. Introduction. In the last chapter, we discussed a special | 
routine Рог "forced first collision." Conceptually, this routine con- 
ducts a particle from a special point of departure, namely, the source, 
to a point of collision within the system, with no escape alternative 
for the trajectory itself. When this (Во) routine is used, it is 
entered only once, namely, directly from the source, and has only one 
exit, the collision routine (7). 

In all problems, whether this device is used or not, a routine 
(B) is required which is designed to conduct a particle from a perfectly 
arbitrary point of departure in an arbitrary zone, at which the particle 
parameters are known, to its next point of collision within the zone, 
or to its point of departure from the zone, in the event the boundary 
is reached without collision. (Cf., however, 89 of this chapter.) 

This (В) routine is entered directly from the source, in prob- 
lems not using the (Во) device, and in various other situations, namely, 


as a particle departs from a collision, as it enters a new zone, or 


-63- 


Google 


upon re-entry into the same (central) zone, after crossing a central 
hole, etc. 

Aside from the initial determination of the distance А to 
collision (calculated for an infinite homogeneous medium), this routine 
is purely geometric, and consists essentially in comparing the distance 
with the distance а to the boundary of the zone along the line of 
flight. If Yl «à, space and direction coordinates are set at the point 
of collision, and one proceeds to (У). If ) 2d the particle is con- 
sidered to reach the zone boundary. In this case its parameters are 
set at the boundary point, and one returns to (В). 

It is clear that when & problem involves zones of different 
geometric shapes, several such routines (8,), .$- 1,2,... may be 
required, each designed for the geometric problem of its own type of 
sector. In such cases the particle carries an additional parameter 
which indicates the type of geometric zone it occupies at any given 
time. Transfer is made from source points, points of collision, and 
So on to the appropriate (8,)- 

One may also note that, under our procedure, & particle reaching 
the boundary of а zone is referred back to (В), unless escape from the 
System is involved, and is then treated anew relative to the zone 
entered. An alternative procedure using & single random number to 
determine the eventual position of collision, or escape from the system, 
by reference to all path segments defined on the line of flight by all 


zone boundaries may be used but seems clumsier to handle, and we do not 


Fie 


Google 


discuss it. The method involved should be clear from the discussion 


at the end of §4, Chapter III. 


2. A routine for the spherical shell. Consider the problem 
of a particle with parameters R, w, E, g, 3. at а point of departure in 


а shell of radii В The flow diagram of Fig. 17 will be seen 


< В.. 
23^» 
to keep computation at & minimum. Here, R, is the radial distance to 
the point of collision in an infinite medium of free path M t is the 
tangential distance to the inner boundary from the point of departure 


at В, and м, = -t/R is the cosine of the angle 7, from OR to this 


t 
tangent (cf. Fig. 18). Reference to the cos у curve of Fig. 19 shows 


2 


that, w being negative, w,^ - w^« O implies that the line of flight 


t 
cuts the inner boundary. Moreover, if Ry © - - 
sign of f E t* distinguishes between & point of collision on one or 


is non-negative, the 


the other side of the inner boundary. 

It will be seen that for "solid spheres" with no central hole 
this flow diagram will work automatically if 2 = ] is the index of the 
innermost zone, provided Ro = О be stored together with the other Zone 
radii R,«...«R,. Thus no special (8) routine is required for the 


central spherical zone. The exits (»') indicate collision within zone 


(cf. the next section). 
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Fig. 17 
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Fig. 19 
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(9) 


3. Reorientation formulas Рог the spherical shell. In case 
of collision (”'), one provides a routine for computing the parameters 
R,w as they exist at the point of collision, before proceeding to (v) 
itself. This is indicated in Fig. 20. Moreover, in the event of escape 
at Ry? we have two alternatives. If у= » ; the number of the outer- 
most zone, the particle leaves the system, and one proceeds to (є), 
while if y», the particle enters а new zone э + l. Classification 
of an escaping particle may or may not involve its direction of escape. 
In the latter case, one may by-pass the reorientation part of the escape 
routine of Fig. 20 by putting the y+ l3 and »- Y boxes first. Опе 
obtains the w' formula of the latter routine from the relation w' 
= cos 7' = ү E - sinf "| and the law of sines: sin "Ry = sin y'/R. 
Note that cos У! takes the positive square root since an entry to a 
zone from the inner boundary has acute angle of entry. 

Similarly, one obtains w' - -ү E - (R/R ^ Q. - «| for the 
entry cosine at Aye For solid sphere problems, the escape contingency 
at R .] occurs only if ү > l, and involves the substitutions w! 3 м, 


using the latter formula for w', H 1 


> R, у - l> vz and thence passage 
to (В). 

If the central zone = і is vacuum, escape at Н, from zone у= 2 
necessitates re-entry of zone у= 2 after crossing the hole. Thus for 


the hollow sphere problem, we have the routine of Fig. 21. 


(9) 


All the geometric relations involved in the present section are 
indicated in Fig. 22 and Fig. 23. 
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Fig. 21 
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4, Flux problems in spherical geometry. Certain problems, опе 


of which is that of a central neutron source in homogeneous air, require 

study of the number of neutrons of energy E which traverse (imaginary) 

spherical surfaces at varying distances from the source. Such situations 
may be handled by subdivision of the medium into spherical zones, and 


т 


incorporating into the "Escape at Ray and Ry-1 routines a cumulative 


tally in counters N These are not terminal categories, except for 


P8 


4-7 and do not enter into a "sum check," but count all neutrons of 


all energy groups whenever they cross a spherical boundary Ry 


5. A routine for the finite cylinder. We consider now the case 


of a particle with parameters x, у, Z, u, у, W, E, g at some point of 


departure within a finite homogeneous cylinder of radius R., height H 


1 
(cf. Fig. 24). The procedure is indicated in Fig. 25. Note that before 
entering (7) in the event of collision, all parameters are stored in the 
machine as they obtain at the point of collision. The direction coor- 
dinates ч, у, у аге the same at the entry to (7) as they were at the (8) 
entry, and do not require new evaluation, as was the case for w in 
spherical geometry. The exit (є) denotes, as usual, escape from the 


system and leads to a classification routine for escaping particles. 


Such routines are discussed in a later chapter. 
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6. The finite cylindrical shell with central hole. Let 


х, у, Z, ч, у, м, Е, g be the particle parameters as they obtain at 
an arbitrary point of departure within a finite cylindrical shell of 


radii В < R and height Н. The procedure of §5 may be used unless 


1 
the line of flight x' = x -ut,y'- y + Vt, 2' = 2 + wt, t20 
cuts the (infinite) inner cylindrical surface x* + уг = n ^ at two 


real,distinct distances 


e 


where $ = ux + vy and 4 = 6° - (1 - м2) (х2 +y - BE In the 


latter case one provides an additional routine as indicated in Fig. 26. 
Study of the flow diagram together with Fig. 2] should make 


the method clear. Note that the case w - l and the case of & line 
of flight tangent to the inner cylinder are handled automatically. 
Nevertheless, some caution -— be required at entry to the "t" boxes 
if l -w° is very small. Such а case may be handled easily by а pre- 
liminery comparison of | -o * Vi, - K(1 - м2) with О, where К is 
a constant larger than the largest distance within the shell. In case 


the latter difference is positive, one may consider EN as essentially 
infinite and may route the flow to the proper box directly, by-passing 


the explicit computation of + and 2... 
An alternative method which is very convenient when several 


cylindrical zones are involved is that of Fig. 27а. This has the 
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Fig. 27 


advantage of effecting the decision between escape or collision with а 
minimum of square roots, sets up the space parameters at escape position 
with little repetition of code, and avoids the "infinity decision" re- 
ferred to above. It is understood,that the inner and outer radii 

о > Bes В, and base plane z-coordinates H' апі Н" are properly set as 

a particle enters one of the system of zones. Moreover, a parameter 


s X + s is carried throughout with x, y, 2, and an additional 


R 


2 


р 1 - у? together with п, у, м. Ап index Х is also used, being 


set equal to zero at the source. The latter is of a purely computa- 


tional character. 
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Fig. 27a 
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It may be helpful to note that q = -1 ога = +1 according to 
whether the directed line of flight does or does not cut the inner (in- 
finite) surface, while d is "plus" or "minus" according to whether the 
point of collision is outside or inside the (infinite) base planes. 
After first transit through the box in which X, 9 Yi» 2, RY are com- 
puted, these quantities refer to the point of intersection of the line 
of flight with the base plane z = H, when d 20, and to the point of 
collision when а < O. The exits Pa and Вр indicate escape from the 
zone through the inner or outer lateral surface, and from one of the 
bases, respectively. In an actual problem involving several sectors, 


such exits must lead to rather involved routines for deciding the sector 


next entered and setting up its geometric parameters RS, Ro 


1? H', Н", and 


yor classifying escape if such is the case. 


Т. The spherical shell in absolute space. Even when the medium 
is spherically symmetric, it may be necessary to keep track of the direc- 


tion of a particle in absolute space, e.g., in case emergent particles 
are to be classified with respect to their directions relative to a 
given source direction. In such Р it is convenient to use х, у, 

2, Uu, у, W parameters. The — is similar to that of the preceding 
section and should require no further explanation. The procedure fol- 
lowing a "core hit" depends on the problem апа may involve an absorption, 
passage to an inner zone, or crossing of а central vacuum. We do not 
elaborate this case further. The procedure is indicated to some extent 


in Fig. 28. 
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8. Slab geometry. In problems on plane slabs, the medium can 


be considered to occupy the region of x, y, z space defined by Z $ 2 © 


= y 
the coordinate 2 and the cosine м of the angle у which the line of 
flight makes with the positive z-axis being the only relevant coordinates. 
The slab may consist of zones of different kinds of media, with upper 


boundaries defined by Zs <... < Z_. The (8) routine Рог such a problem 


may be like that in Fig. 29. 


Fig. 29 
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9. Problems run in cycles of time AT. Іп problems involving 
the distribution of particles (in energy, position, and direction) at a 
given time 4 т from their origin in some source distribution, a parti- 
cle is followed until it is lost to some terminal category such as es- 
cape, capture, etc., or until the allotted time 47 has expired. Such 
problems arise naturally in a-determinations (cf. Chapter II, 88), the 
computation being performed in successive cycles, the output distribu- 
tion of each cycle being used as the input of the next. The (№) routine 
in such cases must be modified so that the essential decision between 
collision within the zone or arrival at the zone boundary is contingent 
upon the qualifying condition: "if time permits." Thus, the distance 
d to collision or boundary must be computed, and the time (cf. Chapter 
II, 83) 

т = т + k"àa/Vk 
of these events is then compared with cycle time AT. If т! < AT one 
substitutes r'—»r and proceeds as usual, whereas, if r' 2 Ar, time 
runs out before the event can occur. Опе then computes the distance 
l = (ат-т)к' УЕ 

that the particle can travel in the time remaining, and the position 
and direction after this distance is traversed. The particle is then 
classified in an energy, position, and direction category "gi and 
one returns to (a). We include in Fig. 30 an example for а homogeneous 
sphere with radius В;. See 552,3 of the present chapter for (у!) and 


7 


Chapter II, 88, for (a). (Lp refers to loss from the boundary.) 
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СНАРТЕВ У 


THE COLLISION ROUTINE FOR NEUTRONS 


1. Introduction. We have seen in the last chapter how the 
geometry of the system determines the immediate fate of a particle, on 
the basis of the equation l= - An r, as a collision or an escape at 
the boundary of the zone. The present chapter is devoted to the methods 
involved in dealing with the former contingency in the case of neutrons 
colliding with nuclei. In the following chapter, we consider the collision 
routines for photons. 

It may be that the media occupying different zones т" 1, 2, TP d 
are so diverse in the types of nuclei contained, and thus in the types of 
neutron processes involved, that it is not worth while to attempt а gen- 
eral code for (7) covering all contingencies. In such cases, different 
collision routines ra) may be provided, each economically adapted to 
its own type of medium. 

The basic object of this chapter and the next is to show how, 
after collision of a neutron or photon in а given medium, the new E, в, 


v, and М are obtained, as well as the cosine of the laboratory angle of 
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deflection from the incident line of flight. The exit (5) from the 
(у) routine refers to a purely geometric procedure which determines 
the new direction parameters м, v, W (or simply w) from the incident 
direction parameters and the laboratory deflection cosine, and thence 


leads back to (В). The (8) routine is developed in Ch. VII. 


2. Capture and selection of the type of collision. It is impos- 


sible to give a perfectly general procedure for the (7) routine, so 
diverse аге the various types of processes. Басһ problem must be studied 
in its own individuality. 

If the number v of collisions is among the neutron parameters, we 


may begin with 


Ф>Б >> 


Unless the medium in question consists of only one type of nucleus, 
and that nucleus has only one type of cross section for neutrons, we 
must next proceed to decide on the type of nucleus hit, and the type of 
collision. Recalling the discussion of cross sections (Ch. 11Т,$$1,2,3), 
it is clear that the total area presented to a beam of neutrons of energy 

| : : А 

group g in the thin slab there defined is (н 2 (tot.) +...) adl. 
Hence, assuming a collision in the medium, the probability that the 


collision be with a nucleus of type A is the ratio of areas: 
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A » (tot.) a a$/ (N А (tot.) +...) аад 


where the «аў may be cancelled. 

If one or more of the nuclei А, B, C, ... present admits 
capture, or some other process which may be regarded &s terminal for 
purposes of the problem (e.g., inelastic collision in monoenergetic 
problems), we may decide to use a weight parameter W, as discussed 
previously. (Cf. Ch. II, 82.) It is then economical to store deter- 
ministic fractions e of weight captured on collision, where 

С = [я А (сар.) + Mn, o? (tot.) + ... | 
the summations being over nuclear types А, B, ... . We should then 
include at the outset the routine of Fig. 31. Here Г. is a "terminal" 
category in the sense that capture is a final event in the life of a 
physical neutron. However, the use of weights prevents the loss of the 
geometric path being followed and greatly improves the statistics. 
When weights are used, it is advisable to use a "weight cutoff" Wo» 
below which weights are negligible, апа an additional terminal category 
L to catch weights falling below the cutoff. This category is ter- 
minal in the mathematical sense that one returns to (а) for a new 


source neutron in case of such a loss, the trajectory terminating at 


this point. If the weight W of the uncaptured beam exceeds Wo? we 
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proceed to decide, if necessary, upon which type of nucleus А, В, С, 
the scattering takes place, probabilities being now dependent on the 
assumption that a non-capture collision occurs. 

suppose for simplicity only two nuclear types A and В are 
present. For machine purposes, we assign to type A a "type parameter" 
value е = 1 and to В the value е = 2. We store in addition to the 


above с, the probabilities 


Fig. 31 


Google 


A, = НА (92 (tot.) - от (cap.))/IN, (o^ 


Р ( tot.) 


= » (cap.)) + Na Gi (tot.) - » ( cap. ) )] 


for non-capture collision on type A, and enter the decision routine of 
Fig. 32 for determination of the parameter е. 

We then proceed to routines designed to decide the type of colli- 
sion on А (е = 1) ог оп В (е = 2). If the processes involved in the 
two cases аге sufficiently similar one may join the exits of Fig. 32 ав 
shown and go to & common routine which, by use of the variable e, can 
handle both cases. Otherwise,one provides separate routines for the two 
types. In any case there is at this point a real disjunction, the neutron 
hitting one or the other type of nucleus. 

It remains to decide which type of collision is undergone, assuming 
it to be non-capture on nuclear type e. This involves reference of a 
random number to an additional set of probabilities, the number of which 
will depend on the number of kinds of processes other than capture which 
nuclear type e admits. 

For example, if e- l and A is uranium, we might have three 
possibilities: elastic scattering, inelastic scattering, i.e., (n-n) 
reaction with loss of energy, and fission. We should have to store in 


this case the probabilities 


e - А (е1.) /le; (tot.) - ^ ( cap.)] 
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А А А А | 

= el.) + o (in. 9 (tot.) - 0 cap. 

1 = [oh Cer.) + of (за) Got - of Ceap.) 

and proceed from the box of Fig. 32 to the flow diagram of Fig. 33. 
It may be preferred, even if capture is to be treated by weights, 

to first decide on the type of nucleus hit, by reference to the total 


probabilities 


A 
№, в tot. 
д бр (tot. E, 


A B | 
LAM (tot.) + LAM (tot. )/E.., etc. 


and then to store individual capture fractions 


г (сар. )/ов (tot.) 


for each type of nucleus, allowing the uncaptured weight to undergo one 
of the remaining processes by use of probabilities such as Es and 3. 
above. 

If capture is not treated by weights, but is regarded as an event 
terminating the tr&jectory, the procedure of this section is modified in 
the obvious мау, and we do Hot include it. 

In any event, &t the present stage of the flow diagram, we should 
have effected the decision of the type of nucleus hit, апа the type of 
collision undergone by the incident neutron, with assignment, if necessary, 


of а nuclear parameter е, апа adjustment of the neutron weight. 
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The possibilities include, of course, return to (@) if capture or weight 
cutoff terminates the career of the neutron in question. 
We proceed to discuss routines for various individual processes 


commonly encountered in practice. 


ка] КУЖ 
O © 


Fig. 33 


3. Elastic collisions in general. We discuss in this section the 
collision between two particles, assuming conservation of momentum and 
energy. We reserve capitals for vectors, and make the following conven- 
tions: R = (x,y,z) denotes the position vector of a point in the lab- 
oratory system, V = В = (х,у,2), the velocity vector of the point R. 


If Ri = (x, 534524) are the position vectors of a set of point masses 
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j #=1,...) п, and m- рш, is their total mss, then В = m mh, 


defines the position of the center of mass of these particles. We also 


recall that if А = (24; T в), В = (b,, Dos ba) are any two vectors, 


177 А = ra, » the normof A is 


defined as |А| = yar, and the cosine of the angle Ө between A and В 


their inner product is АВ = Za 


is given by cos @ = AB/|A|- JBI. 

Now consider two point masses m and n, with laboratory velocity 
vectors Va and Vos respectively. The total momentum of the system is 
the vector P- Ги. У, » its total kinetic energy is the scalar 
k = 25 m, ү; and its total mass is m= Zm. Since the equation 
mV = rm, Vi holds for the velocity V of the center of mass, we have 
&lways P - nV. 


It is customary to define the velocities, total momenta, and kinetic 


energies, relative to the center of mass, by the equations 


t = = 
У. = у. У 

t ._. t 
Poes zn, Vi 

' -pł ‚2 


We find from these definitions the following relations between 


absolute and relative quantities: 
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P= >ш. У. = Еш. (У+\!) = mV + Р! 

dq ^E 1 L 
xs АЕ т rye < Е „г i i 
k= 55 м, 1 = 5 Еш, (У+У) = sm + VP! +k 


or, since Р = mV, 


t ow nes 1 1 
Р' = 0 = my Vi + mo Vo 
k -inv + К! 


Note that the first of these equations states that, relative to 
the center of mass, the particles always travel on the same straight line 
with oppositely directed velocities. 

We assume now a collision between these two particles, adopting 
the notations indicated below for the relevant quantities before and 


after collision: 
Before After 


ү, ГА 
Р Q 
k р 
у W 
ү! W! 

Д: 
Р' Q' 
k' A' 
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and restate the fundamental relations 


Р' = 0 = 8! (1) 


у д==ш\ + д (3) 


Since Р = 9 ала Р = шү, Q= m4, we have У = W, which states that 
the center of mass of the system proceeds unperturbed by the collision. 
This equality, together with the equations (3), and the fact that k = £, 
implies that k' = £', that is to say, the relative kinetic energy is 
unchanged. 

Now Р' = О in general, and we have therefore the HE UH 


2 2 2 2 А " 
m Vi = = mV, nV. = nV » and mV! - m, 2 = О, with identical 


l 


results for the relative velocities Wi, W5 after collision. 
Hence the equations 
2 1 2 
= t = ү Уыз: Е 
олу wq о. ж 
Q. qoc SN 
ту MI - n, Vi = О 


Google 


2 2 
are satisfied by v^, va and also (since k' = £') by Wi and W^ ? 
Since the determinant of the linear system is - Ц ae _ ) #0 

уш ыс Шш с: › 
2 2 2 2 
æ » . 1 Me ' t d 1 
the solution is unique and Wa = MI 7 W5 = V5 | 


Ме have, therefore, the following very simple picture (cf. Fig. 34), 

The upper part of Fig. 34 indicates the relations we have derived: 
The center of mass velocity is unchanged, the relative incoming (and out- 
going) velocities are collinear and oppositely directed, and the relative 
speeds are unchanged. The dashed lines are coplanar, as are the solid 
lines; however, the two planes they determine are not necessarily the 
same. 

The conditions for elastic scattering do not determine the angle 
of scattering, as defined (say) by the angle y; from М to Wi (nor 
the angle between the two planes referred to above). This must be given 
in the form of a distribution law which depends in our case upon the 
neutron energy and the nucleus hit. Such distributions will be discussed 
in the next section. Our immediate object is the derivation of formulas 
for the angles Vi of scattering in the laboratory system, that is, the 
angles between W and the М and for the kinetic energies 
LA = =m, М of the two particles in the laboratory system, as functions 
of an arbitrary angle Vi. 


We first consider the м: 
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we = (W+ ws) = (ү + wi)? = ү? + ZVW) + wie = ү + vie + alv| м | сов "1 
= ү + vj? + e|v||vi| cos v1 
Wa = (W + WC = (у + Wa)? = v^ + omg + wee = v^ ew. +2] | [М | сов фр 
= V^ + vse - elv| KA cos ya 
Next, we compute the cos E 
cos V, = ww. /|v| KA = V(W + му) |у W, | =V(V + 93/ |v | vl 
= (y? + vi! )/|у| [м.| = (vê + |v| [wy] сов »:)/|v] [4] 
= (|У | + |у. | cos v1)/ |м. | 
соз №, = W,/|v||W,| = (м + ма) ум | = v(v + w3/|v] м 


(v? + wi3)/|v||w,| = (v? + Iv|| wa] сов v 2)/|v| |w 


2 


(|v] - Iva | cos v+)/|w,| 


where the СА 2 are given above. These formulas provide the general 
solution to our problem. In most cases the energies of neutrons are 
sufficiently greater than the (thermal) velocities of nuclei to admit 
the assumption that the latter are effectively stationary in the 


laboratory system. 
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We proceed therefore to specialize to the case of neutrons of mass 


ш, laboratory velocity V. scattering on nuclei of mass n, velocity 


1 
V5 = О. In this case we have 
_ " | -l 
Ур diae ec 43 MEE m, |v, | 
Vey -vest oam V. ea aV 
j 1 “| 2'1 
| = ENS. NEZ 
V5 = У -V = V = -m mV. 


Substitution gives us 


2 22 


W = 2 (27 ^u^ + ш шр + 2m ^ Im, — 


. Q1 2 T 
and hence, if f = 5 CENA у К, = 5 m V. are the (laboratory) energies 


of m after and before collision, we obtain 


А/к, AE: + 2m m, + ms - emm, + emm, COS vs} 


1 - 2n ^nm, + 2n ^nm, сов y: 


It is convenient to introduce the quantities A = ш o/m , 


(A - 1)9/(A + 1)* "" - qm, )*/(m, + п, jeter ^n m, 
S 


1 - Bru -2 - 
5 (l+r)=1- am ^mm, and Т = (1 г) = 2м. mm,, obtaining 


finally 


£ - ki (S 4T cos фу) 
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Moreover, 


-1 ai | 
cos V, = [v {п m, + ш m, cos Vi Ум, 
= (575, * mm, cos 7 ү {22 + m “ш + 2n ^nm, сов vf 
\/ 2 
= [1 + cos 737 fsa + 2А cos 91) 


Although we make no use of them in the sequel, we include for 


completeness the results for the scattered nucleus m 


vI 


ХЕЛЕ = zc - cos и} = 2T sinf (3) 


and cos P = ү 5 (1 - cos Vi) - sin (1) ; Which the reader may verify 


as an exercise. 


The salient results of this section are the formulas 


E'-E(S,-T л) 
e e 


а, - {1 + AL H ү + Gs + Аи} 


where we agree on the following definitions: 


m mass of nucleus of type e 


m= mss of neutron 


Ao = m /m 

rou ee 1)*/(&,. + 1)° 
Б = 3 (1 + г.) 

Te «$ ео 


Е = neutron energy before scattering 
laboratory system 
Е' = neutron energy after scattering 
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И = cos yi 


Vi - angle of deflection in center of mass system from original 


line of flight (coincident with V - W since узе О) 


а = cos ul 


ф = angle of deflection in laboratory system from original line 


of flight 


The constants Ao? 5е Те are nuclear (energy independent) constants, 
the subscript e, if required, being set in the routine of the preceding 
section. 

It may be noted that in the case of scattering on hydrogen, we 


may assume А = 1, and the formulas become simply г = 0, S = Т = 5, 


Е! = 5 E (1 + и), anda = ү5 а + и). Several remarks are of interest 
here: (a) we have the relation E' = Ea‘; (Ъ) scattering on hydrogen is 


always forward in the laboratory system, viz., a = cos i 2 0; (c) 


а = COS 9» - ү х (1 + cos $1) = COS (v. /2) implies Д = 41/2; (а) е 
formula а = 5 (1 + #) should be used for hydrogen, since the general 
formula is indeterminate at u= -1; (e) if scattering on H is isotropic 
in the center of mass system (a good assumption, incidentally), we have 
#= 2r - l, and a= yr, which shows that the laboratory scattering is 
in the cosine distribution. (Cf. Chapter II, 850.) 

Finally we observe that for "heavy" elements (4, large) we HEN 


the following approximate results: f. 7-1,8,-21,T7,-0,E'-E, a - n. 


e 
If these relations are of acceptable accuracy, one may ignore the energy 


change and determine a directly instead of M by the methods of the 
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next section. In this connection one should refer to Chapter VII, 65, 


for the case of scattering isotropic in the laboratory system. 


4, Тре differential elastic scattering cross section. Ме have 
seen in the last section that the preservation of total momentum and 
kinetic energy characterizing the elastic scattering process does not 
serve to determine the angle vi of scattering in the center of mass 
system, but does determine the new energy Е' о? the scattered neutron 
and its deflection cosine a in the laboratory system in terms of a 
given И = сов Va 

The determination of И is governed by a physical distribution 
function e; (9), whose units аге barns per steradian, and which depends 
upon the scattering nucleus e апа the energy Е of the incident neutron. 
Specifically, T. (9) à9 is defined as the cross section (in the sense 
of Chapter III, 81) presented to incident neutrons of energy E by a 
nucleus of type e for the process of scattering in the center of mass 
system at an angle in the 49 neighborhood of the direction 2 with the 


incident line of flight. Thus by definition, 


IE: (9) a =o (el) 


the integration being over the entire surface of the unit sphere in 
direction space U, V, W, and oF (е1) being the (total) elastic scatter- 
ing cross section for the element e at energy E. In terms of spherical 


coordinates (azimuthal angle g, polar angle Vis Fig. 35), we have 


Google 


Google 


incident | line of flight 


Fig. 35 
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deflected line in C. М, system 


т 2х 
ое (9:,Q) sin v; афау! =ое (е1.) 
-— SUR 1 1 Е 


In all problems of the text о = 


p 18 independent of @ and we 


have 
2л | о (41) sin Vai d Vi - ор (е1.) 


or, in the form usually given, 


| 


L 
e e | 
erf e" (и) ам ‘р (е1.) 


It follows that рр (4) = 2л ор (и)/ ое (е1.) is the probability 
| 
density function for elastic scattering at the direction di from the 
line of flight in the center of mass system. 


Hence, the Monte Carlo procedure sets 
u 
r= | M (и) du (-1 3 u € 1) 


for the determination of и from the random number г. 


-101- 


Google 


5. А routine for elastic scattering. Suppose that elastic 
scattering on the (light) element e is isotropic in the center of 
mass system. This means that the function c^ (и) is а constant, and 
the final formula of the preceding section gives simply #= 2r - l. 

In such a case, we should enter a routine like that in Fig. 36. Here 
Е. represents the "energy cutoff," namely, the last of the lower bounds 


G 


E, > E, >... OF 


1 of the energy groups adopted. Energies falling below 


G 
this necessitate classification of the corresponding weight in a terminal 
category Le reserved for loss to energy cutoff. The loop on = - Е 
begins by comparing the new energy Е with the lower bound Е of the 
group which the neutron occupied before scattering, since elastic 
scattering cannot raise the energy. Note also that E, - Е < O at entry 
to this loop ensures its termination at some g 2 G. If Ес > О, the 


previous E - E, decision is mandatory. The formulas for E' апа а 


G 
are fully discussed in the preceding 83. Аз remarked there, the pro- 
cedure for the special case of hydrogen is simpler and is indicated in 
the lower half of Fig. 36. For the case of scattering isotropic in the 
laboratory system, we again refer to Chapter VII, $5. This applies 
approximately to the case of heavy elements with an isotropic law in 
the center of mass system. 

The (4) of Fig. 36 refers to & routine for fixing the final 


laboratory direction parameters of the scattered neutron (cf. Chapter 


VII) before returning to (f). 
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(е not indicating — 
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Most elements have more complicated differential cross section 
functions, usually characterized by preferential forward (and sometimes 
backward) scattering, an effect which tends to increase with higher 
incident energies. The problem of choosing И from such a distribution 
is essentially that (already discussed in Chapter II, §5) of selecting 
the direction parameter w froma given source distribution. 


If the problem can be restricted to a single energy, or to 


(10) 


relatively few energy groups over each of which the differential 


function 0. (и) may be considered constant, one may well store tables 


for the probabilities of scattering at cosines и 2 H., where 


P 
gJ 


M >И >... ти. = -l are suitably chosen bounds of subintervals of the 


l 2 J 


cosine range, and the 


are pre-computed by numerical integration and stored. The routine for 


determining 4 is then that of Fig. 13, where we read j for i, Pa j 
2 


for P,, А, for w,, and А for м. 
i д 1 


Aside from its demands on storage space, this is subject to the 


errors of interpolation, which may be difficult to minimize in cases of 


(10) ese need not coincide with nor be as numerous as the energy groups 
reserved for free path and scattering-type probabilities. It is not 
unusual to carry two or even three different sets of energy classifications 
with corresponding indices g, h,... аз neutron parameters. 


-l10h- 


Google 


very forward scattering. 
If the original curves оъ (и) can be simply fitted, the simplest 
scheme may be the von Neumann device (Chapter I, §5). We should then 


have a formula for the function 


HA 


) 


% 
о (е, Е И), 


me, В, и) = ор G/ { пвх ор (и) on 1 Зи 


and use the routine of Fig. 5 with а = -1, b= 1, p 


reading и for $. 
* 
Occasionally, the function 9 (e, E, и) may be fitted easily for 


each energy Е ава function of И, say, 


* 
о (e, E, и) =A, + Bou + се ue -1<н<]1 


but the coefficients may prove intractable as functions of Е. In such 
а case we may be able to store the coefficients 


e е е 
А”, В”, С 
& ge? 8 


for each of a reasonable number of energy groups, апа use the routine 
* 

of Fig. 5 as indicated above, the machine computing о (е, в, и) = 

AS + BA + C," 2 from И and the stored coefficients. 


Finally, and in the worst cases, it may be necessary to store 


х 
а table of values of © = o  /(maxo , Where c = 
m,n myn’ ии П m,n 
2 


ALT E) is the differential cross section for & particular element 
(we drop the index e temporarily) evaluated at а suitable set of 


values н = Ioa. а? = -1 апа Е > 2 Eg where E is at least 
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the highest energy encountered їп the problem, and ЕҢ = Eo the least 


such energy. The energies and cosines used for this purpose may be 


handepicked, and need not coincide with those used for other purposes 


elsewhere in the problem. We store a table of the form 


Fig. 36a 


and resort to the double interpolation routine of Fig. 37. 
However we may determine А , the routine should lead to a 
determination of E and g after collision and the deflection cosine 


а for the laboratory system as indicated in Fig. 36. 
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Fig. 37 
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6. Differential elastic cross section for the laboratory system. 
If a differential cross section e. (a) is given for elastic scattering 
relative to the laboratory cosine а = cos Vi one may choose between 
the following two alternatives. 

(а) One may employ e. (a) directly to determine а, just as 
the center of mass differential cross section с(И) was used to determine 


и. It is then necessary to use the formula 


2 
"ES Ca - a^) + (sgn a) V (1 - а^) - (1 - а?) + а2д? } (4) 


to determine the corresponding center of mass cosine и, which is needed 
in computing the new energy Е' = E(S + Ти). 


The above formula is obtained by solving the equation 


a = (1 + Ам) / у [1 + A? + 2ди] (5) 


for И in terms of а. The (sgn a) choice of sign on the radical is 
dictated by the following considerations: From equation (5) it is clear 
that the sign of а is always that of 1+Au. In solving equation (5) 


for BM we obtain 


jd AU cR d Vita - в2)^ _ (1 - в2) + а2д2] 


зо that the sign of the right side must always be that of а. It is 
clear graphically that for the function f(a) - а^ +t V g(a) to have the 
sign of a, it is necessary and sufficient that one use the upper branch 


Рог а > О and the lower for а < 0. 
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Formula (4) may be rewritten in the form 


2 / 2 
~ +а|/1 - 25 

A 
where b? =] - а^ " 


(b) Alternatively, we may use the defining relation 
olt) аи = оү (а) ав, 


to compute the center of mass function о (и) in advance, and then use 


the usual procedures of $5 directly. This may be done by means of the 


formulas 

о(и) = « (a (#)) са 
Mene а(и) = (1 + Au)/ V [1 + А? + 2an] 
and 5а =[А (А «3]/ [1 + А? + оди] 3/2 


T. A weight. device for elastic scattering. In problems involving 
many collisions in a medium admitting elastic scattering on light elements, 
it may prove worth while to avoid loss of trajectories, which are usually 
followed through many collisions, to the energy cutoff Eo. These losses 
may be obviated by the following device, which makes further strategic 


use of weights. 


We have seen that the relation between emergent and incident energies 
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in such collisions is E' = Е(5 + Ти), where и is the center of mass 
deflection cosine, chosen randomly by the equation 
1 
r= | p, (4) ён = B (и) 
u 
Now, for a given incoming energy E > Ea? the least emergent energy 
is E'n” E(S - T). If this exceeds E,, the energy cutoff for the 


problem, one proceeds &s usual, with no possibility of loss to energy. 


ni 
t 


However, if E(S - T) < E,, that is to say, if Е < E/(8 - T) 
< | < 
there exists a critical cosine He? depending on Е, such that 


E(S + Tu.) = E,, with -l< 4. <1, and 


де 
f = f(u.) -f pg (4) de 
-l 
4 
represents the probability of scattering with 4 on -l< а < и with 


the resulting energy E'< E,. Moreover, assuming и > И 


С с 


1 l 
"| pg (и) аи/ | pp (и) аи = 9р (и) 


Mc 


is the proper Monte Carlo formula for determination of the cosine 
on the range и < и < 1, with a resulting energy Е! 2 E We may 
therefore proceed as in Fig. 38. 

For the simple case of scattering isotropic in the center of 
mass system, we have the formulas f = 5 (1 + и), и = qz (r) = l-r(i-#), 
uz Po (г) = 1 - 2r. * 


Finally, we note that & neutron being followed by this method a 
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automatically has weight М > О and energy E OE if the random 


G? 


number г = 1 does not actually occur in our sequence. 
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8. Fission. Suppose that a neutron causes fission of some 
(heavy) nucleus, resulting in the production of an average number A 
of neutrons with energies of release independently distributed in a 
given fission spectrum, and whose emergent direction distributions are 
isotropic in the laboratory system. 

It is possible to invent a probability distribution p(n) with 
= Xn p(n), decide on the number n of progeny in а given fission 
by chance, and follow these n neutrons individually. This may be 
done even in supercritical systems if recourse is had to time cycles 
and census taking. However, this is unnecessary and undesirable. We 
may instead follow & single neutron, of weight £^ times that of the 
incident neutron, chosen from the fission spectrum of energies, and 
directed isotropically. This may reduce fluctuations, enormously 
simplifies the code, and decreases machine running time. 

Suppose that f(E) dE is the probability of fission energy E 
between E and E + dE. We choose a set of energy intervals with 


—— 


bounds E 2 Е >...> Ер, where E is the highest significant fission 


energy and E. the energy cutoff for the problem. These E, need 


not agree with bounds of energy groups used for other purposes. Define 
О 
Е =|. Р(Е)аЕ 


We may then refer to Fig. 39. Note that if the E are differ- 


ent from the E used for (say) storing free paths A one must deter- 


mine the index g, beginning the loop with g = 1 since fission can 
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raise energy above the incident energy Е in group g. Тһе exit (6) 
refers to the routine of Chapter VII, $5, for Scattering isotropic in the 
laboratory system. 

In some problems, designed to find the distribution, s&y, in 
а number of cylindrical shell zones defined by radii В. < R,<... < Ra 


1 2 
and heights Ну < H, а Hy? of fissions resulting from a given spatial 
distribution of fissions, fission is a terminal event, and counters 
Nij are reserved for the number of neutrons terminating in fission in 
radial zone i,height zone j. In the event of fission we should then 
follow the simple classification routine of Fig. 40, rather than that 


of Fig. 39. 


Fig. 39 
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Fig. 40 


9. Inelastic (n-n) collisions in general. Using the notation 


and procedure of §3 of this chapter, we now discuss a type of collision 


in which an incident neutron of kinetic energy ky = 5 n V^ = Е strikes 
а nucleus of mass m,, imparting to the latter a known energy of excita- 


2 
tion e, and emerging with an energy L. = à ши. < = Е'. We obtain again 


the fundamental relations 


Р! = (6) 
(7) 


mv + P (8) 


= 


2 


O=Q! 
P = mV 9 = 
nV + k' (= 


м 


The momentum and energy conservation laws now read 
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P=Q 
EM 
however, and we find V - W ав before, but now 
k'= xe 


We find again that the equations 


determine wre and wee uniquely, the difference being only that A' 
is not simply k' but k' - «e. Thus the speeds of departure of ту 
апа I, in the center of mass system are uniquely determined, although 
they are not simply the incident relative speeds as before. The angle 
va must be determined according to some physically determined distri- 
bution law, and from this the laboratory energy and deflection cosine 


may be computed. 


We will derive the relevant formulas. First, solving the above 


system for Ws Wi, Suppose that we first rewrite the system for 
T y2, 
Vi ; ул : 
| 
5 шу! +5 шул” =k' = f' +e 
2 242. | 
ШУ -шруу = 
in the form 


Google 


1 ‚2 1 ‚2 f 
amU;] - 6) + 5 mV - 65) = 1 


2, 2 2, 2 i 
mV; -6)-m(Y5 - є„) = 0 


l 1 
= € € = Е 2s — = 
where e, = 2m, /mm, > = m /mm, and 5 ше, + mé = 6, 
ще) + ms 62 = О. We see at once that 


2 ‚2 "c c 
1 jo! FONS -% 


We deal only with the neutron т). Computing we as before 


in terms of cos va » we obtain 


rye = 92 = 2 + ie 
(и ew) = (у +М])” = У + Му” + ам [\1| cos vy 


2 ‚2 V ‚2 
үче = e 217] (VIN - 6) cos 41 


and for cos Va ; the laboratory deflection cosine, 


2 
"1 


cos фу = W4/|V||W,| = vw + М) у |м = ve + м) УМ | 


1 


(v^ + ун.) / [v] м = (y? + MIEA сов 44 )/|У| |w] 


(|у | yii- cosy )/|W. | 


i ors _ Li. dedi 7l 
Again specializing to the case Va =0,V=m mV., Vi -m mV. , 


(|v] + LM cosy )/ [М | 


we obtain 


we ЗА + ш nV - є + 2n ^n. | V,| cos V1 у (n “neve - €1) 
М + m n (D mtm, соз pj Ya -P Aer) 
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and therefore 
A/c, = 72 {пу + amma + пр (30е AT - amm + асови 1-08 va} 
= (1 - 2m, m, /n' ) = (e, We) + 2n mn ^ COS и! E - (шее, ny? ) 


E 2 = _ 
But € = 2n, e /mm, so that m є fna = шешу = me/m E, and 


hence, finally, replacing А by E' and k, by Е 


1 


Rak {s (3) (me/m E) + T cos vi Vi - ела) } 


Moreover, 


сов V1 = аА + cos vo (m 225ү e M [V,| 


= |v CEN E nm, SOR yy Mie oe (n^e,/ АД! ДШ 


nm, + mm, cos V1 ү: - пек) 
Finally 
1 +A cos pj Va - Gen) 


Now the reaction can occur only if 4 =k' -€ 20; since Е = 


Or 


cos 9i = 


cos Vi = 


1 e co 2 -2 2 -21,2 _ - | 
L3 m Vi = {а + mnm у = mm lg ‚ the condition is simply 
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Е >шш є. If we define this critical energy as є„› we may write 


[=] 
{1 
tj 
—— 
л 
— € 
+ 
"mr 
"> 
w 
№ 
Е ^ 
О 
nasi 
+ 
ЕЗ 
Q 
О 
[42 
€ 
[E 
кз 
| 
m 
Ош. 
Dr] 
— “=” 


1 + cos v a - (5/8) 


cos Vi = Е | ——- 


fı „д^ [1.2 («,/£)] + 2А cos „|: - Cem} 1/2 


About these formulas we make the following remarks: (a) For 
J^ = О, they reduce, as they should, to those Eos elastic scattering, 
(b) for E. > Oand A large (heavy nuclei), they become approximately 
Е' = Е -є and cos Vi = COS Vi as we should expect. 

It must also be pointed out that а nucleus may have various 
excited states each with its own excitation energy є , and а correspond- 
ing inelastic cross section g,(e) for neutrons of energy Е > (m/m, ) е = 
For light nuclei these must be dealt with individually, using the above 
formulas with the appropriate є. 

For heavy nuclei, the states may be well separated in some cases, 
and may be dealt with individually, using the simple relations of remark 


(b) above. If, however, the excited states are very close together, 


one resorts to such methods as those of the following two sections. 


10. Inelastic (n-n) collisions on heavy nuclei. If the states 
of the nucleus are closely packed on the energy scale, it is customary 


to give a table of experimentally determined probabilities for 


Pa h' 
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ап incident neutron of energy group h to produce in (п-п) reaction a 
neutron of energy E' Е, where h' 2 h, and Е > аль > Eq is a suitable 


set of lower bounds for energy. One may then proceed to determine the 
energy E' in the usual way by use of a random number and linear inter- 
polation, as in Fig. 41, which is drawn for the case of isotropic 
scattering in the laboratory system, and therefore exits to (ô) (Chap- 
ter VII, 85). Note that the Ph oh! form a triangular matrix, with P 


non-zero only for h' 2 h, and P 


h,h' 


h,H = 1, һ = 1, ° ө ө 9 Н. 


Some inaccuracy is unavoidable in such a method, since, strictly 
speaking, Р h! is a function not only of the group h, but of the energy 
4*2 


E in this group. 


ll. Inelastic (n-n) collision with Maxwell distribution. Con- 


sider an (n-n) collision of а neutron of energy Е < E with а heavy 
nucleus for which it may be assumed that the emergent energy E' 


probability density function is proportional to 


n(E') = Е' exp(-E'/T) 0<Е' < Е 


where Т = 4|(Е/Е), а being a given constant. 


-119- 


Google 


Fig. 41 


Е 
If A(E) -f E' exp(-E'/T) дЕ' 
О 


we have 


p_(E') аЕ' = (1/А(Е)) E' exp(-E'/T) аЕ' 


pl 


and the Monte Carlo principle involves determination of E' 


and r by means of 


r “| PQ(E') aE' 
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from E 


The solution for E' is difficult апа we use instead the 


device of Chapter L$5. The latter requires the function 


р (E') = n(E')/max N(E') 


the maximum being for the range O< E' < E. Now the maximum of n(E') 
is at E' = T, and we may use this provided T = a Y(E/E) < Е, which 
means that the incident energy E shall exceed a TE. Let us suppose 


for simplicity that e^/E < E., the energy cutoff (as is usually the 


ЕС 
case). We have then 


pŒ) = (È) Bre E/T <1 


for all incident energies Е оп the range Е. < Е < Б and we may pro- 


ceed in the usual way, with Е' = rE, etc. 


12. А combined transfer matrix for fissionable nuclei. If the 
nature of a problem involving а fissionable nucleus does not demand 
keeping the individual types of collision separate (for example, for 
purposes of tabulating capture, fission, etc.) and if the neutrons 
emerging from each type of collision are all emitted isotropically in 
the laboratory saver. a much simpler method can be used to great 
advantage, in place of the several individual procedures indicated in 
previous sections. Let us adopt the following notations, in addition 


to the usual a (cap. ), a (el. ), a (in. ), o (fiss.), a, (tot. ): 
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p(h,h')= probability of scattering from group h to group h' 


(inelastic collision) 


у: = average number of emergent neutrons рег fission 

f(h') = probability of fission neutron energy in group h' 

б ph F Kronecker ð function, with value 1 for h=h', O otherwise 
2 


Now, assuming a collision of a neutron of group h with the 


nucleus, it is clear that the expected number у nt Of neutrons emerging 
2 


from the collision in group h' is 


Ere = {2 (сар.) ожо (61), pi + оу(їл.)р(һ,һ') + 


e, (fiss. ) f(h' y, be, (вов. ) 


Thus, the total expected number of neutrons per collision is 
v= р | Р 7 {o,(er.) + e, (in. ) + o (fiss. ) v. } /o, (tot. ) 


Hence, the probability of an emergent neutron from such a colli- 


sion being in group h' is 


dn,h' ^ ZU 


We may therefore store & table of ‚= 
Á Sh yh Ё qk 
h S Шо В 4" 9 Deeg Hy and но = 0, h= 1,..., Н, as shown 


below,and proceed according to Fig. 12. 


-122- 


Google 


dii cu = " ES Qn, ht-1 7 Sh ne 


Fig. 42 
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13. Collisions shattering a nucleus. We consider in this 
section a type of collision in which a neutron of rest-mass m and 
= О, 


velocity У. enters a nucleus of rest-mass m, and velocity V 


l e 2 
resulting in а shattering of the neutron + nucleus system into а set 
of fragments of masses п, апа velocities W д = 1,2,... , the 
total rest mass n - 2л, of the fragments being greater than m - mtm, 
We treat the mechanics of the collision non-relativistically, using the 
formula 3 пу? Ғог kinetic energy of a particle of rest mass m, and 
velocity V, and we suppose that ш = п at various points of the argu- 
ment. The actual mass-difference п - ш corresponds to an energy 
Е = (n - nen where c is the velocity of light, and in order for the 
process to occur, part of the kinetic energy of the neutron must be 
used to supply this energy €. 

Proceeding as in §§3 and 9 of the present chapter, we have the 


general equations 


(note the assumption m - n), and the conservation laws now read 


Р = 9 
к = fre 
As usual, V = W, апа so К! = уче . Just аз in 89, we find 


that in order for the reaction to take place, we must have k'= mn ^E 2 є 
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or 


E 2 шие 


Assuming such an incident energy E, the conditions on the 


relative velocities мо the fragments аге 


Хам! = f = ng E -е 20 


Duw'=0 
Jd 


In contrast to the situation in§9, these relations do not deter- 
mine the н; if more than two fragments exist (including the original 
neutron). We propose instead to single out an arbitrary one of the 
fragments ny and discover the maximum relative energy with which it 
can emerge. 

We adopt the convention that 5 refers to summation on all 


j J. Thus we define n' = S'n, to be the mass of the remaining 


J 
fragments, and n'W, = S'n W3 defines the velocity W, of their 
own center of mass. Now Wi = W, - W denotes as usual the velocity of 


R R 
the residual center of mass R relative to that of the whole system. 


: е _ 22 ' E ' = t [ 
Moreover, we define м} = W, ША (w; + W) (We + W) м, We for 


all j#dJ. Now the above system becomes 


(9) 


ам: + Di ам, = 0 
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But observe that the following relations hold for the residual 


system: 
Li "= ' ' " = тт! 1 n 
У Daw 2 g V +} 3 n'We de zd (10) 
x: ie ui E ' nye _ 1 е 1 t " id ne 
E'anHi = D's a, (We + Wi) ерам + Wp ам, + у' п, (11) 
Now we know by definition that 
'n W., = n'W 
EE R 
so that 
' n, (W+ Wi!) =n' (М + ММ 
Qin, | 5) ( д) 
and hence 
'n W! = n'W/ 
D nyja oe 
Reference to (10) yields 
пени = О 
2 J J 
and this result with (ll)gives us the familiar 
t 1 je = 1 ' je ' и ne 
Lig awe -anW + Yi 5 aw (13) 


Now to return to the system (9), its second relation together 


with (12)tells us that 
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so that 3 nw, = (ny /n' )(5 пт” ). We may therefore rewrite (13)as 


(1 = t p UM 
у'н = (nj/n' (5 пи") + зау 


Тһе latter, substituted into the first member of the system (9), 


gives us finally 
2 1 
m — t ! — ! m 1 
(m/n') (5 nw pe OS 5 п М" f 


It follows that the maximum relative energy of пу must be 
(n'/n)E' = (n'/m) (mm Е ые) 

Thus, in summary, if the incident energy E exceeds the miniml 
energy mm," < required for the process, the maximal relative speed 
of any fragment J is given by 


Suppose that we now construct the laboratory velocity Ws and 


the corresponding laboratory deflection angle d by the usual addition 


м = ү + W3 = ү + Ws . Reference to Fig. 43 makes it clear that the 


possible end points of the vectors М. occupy the sphere of radius w 


J J 
about the end point of V. Moreover, it is apparent that three essentially 


different cases arise, depending on the relative magnitudes of Wy and 


lvl, or, equivalently (using the above formula for Wy), the relative 


magnitudes of E ande = n'e/(m, -n If we are interested only in 


д): 


neutron masses n_, the latter reduces to т 5€/ (n, - n). 


d^ 
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Сазе 1. WIS | v | or Е< є. The angle s is limited to the 
range O€ T < arc sin (м ЛУ), апа the possible speeds ГА соггев- 
ponding to any given vs on this range are bounded from zero. 

Case 2. Ws =|v| or Е = Є, The angle V, is limited to the 
range O< i 7/2 but for a given such d the possible Ln range 
from O to & maximum. 

Case 3. мр2 |У| or E>e. The angle d has its full range 
о < V < п, and LE for a given р. ranges from О to its upper bound. 

In any case, the maximum speed w, associated with a given 


J 
possible VI is given by the greater root of the quadratic equation 


Fig. 43 


—2 2 e : А = E 
um - ө Е 
Wy Wy + У 2 ГАА cos v Substituting for Wy and writing J 


for the maximal laboratory energy 5 nj of the gm fragment associated 
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with the laboratary deflection angle у д? we obtain 
-1 1/2 1/2 -1 
€ z € t - > = | 
7 [2m (mn E) cos V, | т tm [n'é (m, n_)E ] O 
If we apply this to any neutron пу, we have ny = m, and 


п' = p » 80 that finally 
€ - (2n * VE cos y.) eY2, wll € - (m, - JE ] = 0 
ш “1 J' 4 mo Bom 


The discussion of this section thus serves to show why the in- 
elastic cross section o, (in. ) for such & process becomes zero for 
incident energies E below nme, and how the energy range of 
neutrons emerging from the reaction depends on the incident energy 


and the cosine of the laboratory angle eL 


14. The (n-2n) reaction in deuterium. Ав an example of the 


preceding theory, we consider the reaction 


2 1 
п + уН —» Н + en 


in which а neutron disrupts the deuterium nucleus into its constituent 
proton and neutron. 


Adopting the atomic masses of the following table 


neutron 1.00893 
H 1.008123 
H 2.014708 


we find an excess of .002345. This gives an actual mass of 


Google 


.002345/A gm per deuterium atom. Multiplying by ES and dividing 

by the number 1.60203 x 10-6 of erg per Mev (Chapter II, §3) gives 

є = 2.184 Mev. Hence m ше = 1.5€ is the minimal energy for the 
process to occur, and the three cases mentioned in the preceding section 
depend on the relative sizes of E ang mé / (m, - m) =2е. 


Ме will suppose that а "cross section" T(E ‚ V) is defined in 


such a way that 
2л T, (E', V) dE' ü(cos р) хада в 


is the expected number of neutrons of energies between Е! and Е! + dE’, 
and directions between ¥ and ¥+ аф resulting from such inelastic 
collisions of a beam of B neutrons of energy E traversing a medium 

of numerical density М and thickness АЙ. (Contrast Chapter III, $1.) 


Now 
ea. ( in.) № 4.В 


із the total expected number of neutrons in this traversal, since each 
such collision produces two neutrons. 


Hence we have 


1 p€,(E,y 
2” | | г i TUS y) дЕ! а(сов V) = 2e, (in.) 


where the € ,U ‚ V) &re the upper and lower energy bounds referred to 


in the preceding section. 
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Moreover, 


-— EE (E', V) аЕ' d(cos y 


E 


is the probability of an emergent neutron being in the range (E, E + dE) 
and (¥,y + dV), assuming such a collision occurs. 

The Monte Carlo procedure is therefore clear. We decide in the 
usual way whether the collision is of this kind by reference of a random 
number to (10. )/o,(tot.). If г is less than this ratio, inelastic 
collision occurs. We then double the weight W of the incident neutron 
and determine its cos¥ =a and its energy Е! from the above probability 
distribution. We may do this in two steps, using first the probability 


density function 
E “2 
Pala) = тор Gaf, ТЕ (E', V) дЕ! 
l 


to determine the а = cos №, and then finding the energy E' using the 


density function for E': 
€ 
| 2 
gg (&,E') = 087907], т(Е', V) aE' 
l 


for the а = cosy determined. 
A more complicated example involving ап energy cutoff and 
critical angles in the same way as indicated in $1 of this chapter may 


be found in LA-1606 (not available). 
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15. Ап (n-2n) reaction оп heavy nuclei. As a further illustra- 


tion of the use of weights, consider a collision of a 14 Mev neutron with 
а heavy nucleus which results in the emission of a pair of neutrons, 
isotropically distributed in the laboratory system, one neutron being in 


the energy distribution 


E! exp(-E'/T, ) дЕ, o<eE's< 


IA 
wm 
а 
М 
№ 


апа the other in a similar distribution 


WA 


E! exp( -E'/T,) дЕ, OSE'S є, Tce 


where €, Ti» То are constants. 


Let А, - | в'еко(-в'/т,) dE' = T fı - (1 + т.) exp(- e /т,)} 


for 1 = 1,2. We have then 


p, (E') аЕ' = АТ” Е' exp(-E'/T, ) dE! 


for the probability of the ho neutron emerging between E' and 

Е' + dE'. Since the expected number of neutrons in the latter range, 
per collision, is p, (E' )aE' + р„(Е')дЕ', we may properly assume опе 
neutron of weight W suffering such a collision gives rise to a neutron 


of weight W, chosen in the energy distribution 


р(Е')аЕ' = = [p,(E') + р(Е')] ag: 
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Since the maximum max p(E') on the range Of E'< є can be 
computed in advance, we may use the von Neumann device of Chapter I, §5, 
followed by an induction loop for determining the new energy group, and 


the procedure of Chapter VII, §5,for the final direction parameters. 


16. Capture in a small zone. Certain problems involve determina- 
tion of absorption in а small zone of material, surrounded by a relatively 
large system of moderator, which slows down neutrons by scattering on light 
nuclei. The slim chance of a neutron hitting this small capture zone 
may be improved by various devices, of which we indicate only one. 

Suppose for simplicity that a spherically symmetric system con- 

E Let R > Ri 


If & neutron undergoes collision 


tains & small central core of capture material of radius R 
be chosen comparatively close to Ку. 
at a point of radius R 2 в" ме пау proceed as usual. However, if 

В < в“, we may process from this point to termination m neutrons each 
of weight W/m, instead of the customary single neutron of weight М. 


EN + 
If the multiplicity m is sufficiently large and В - В. sufficiently 


1. 
small, some of the m descendants are very likely to hit the absorbing 
core. 

Such & scheme requires modification of the over-all flow diagram. 
We call &ttention especially to the following changes: (а) the source 
routine (о) sets а new parameter m = О for each fresh neutron leaving 


the source; (b) the (В) or (B) routine exits, in event of collision, to 


а new entry (У) which is preliminary to the usual capture routine (7); 
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(с) the (7) routine, in case the multiplication trick is indicated, 
stores all parameters В, м, Е, g, W,... of the "parent" neutron at 
the point of collision in new vositisnà designated by R, W, E, g, W,... 
for reference in processing each of the m progeny from this point to 
termination; (d) the (а) entry, to which one returns on termination of 
any particle, is modified to order complete processing of m progeny 
before starting out a new source neutron. The parameter m = 1,2,... m 
indicates the number of the descendant being processed. The following 
flow diagram (Fig. 44) includes the essential modifications. 

Note that the entry (>) excludes using the multiplication trick 
on a descendant. That is, we do not iterate the process. In practice 
the device is usually of a more elaborate nature, using different 
multiplicities m for different critical radii в". Moreover, the 
capture cross section of the core usually becomes significant only at 
low energies, so that the multiplication decision may also rest upon 


* * 
а decision on E-E where E is some stipulated low energy. 
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17. Capture by a "point" detector. The method of the preceding 
section may be ineffective if the capture zone is extremely small. The 
present section presents a much simplified example of such a case. 

Consider a spherical shell (radii R} < R5) of a single light 
element which can be considered (perhaps by use of transport cross 
sections) to scatter neutrons isotropically in the laboratory system. 
Suppose further that a detector of very small radius Ro << Ry is located 
at the center of the hole (vacuum), and it is desired to find the dis- 
tribution into energy groups of neutrons impinging on the detector, con- 
sidered as a perfect absorber. 
denote the bounds of the 


Specifically, let Eo >Е.>...> Е 


1 С 


energy groups adopted, Ес being the "cutoff" energy, below which neu- 
trons are lost to a category Lp. Consider a collision С at x, у, Z, 


В, В. & R«R, of a neutron with incident direction м, v, м. The 


1 2 
fraction f of total solid angle subtended at С by the detector is 
2 
= (1 - сово) = = (1 - (1 - if) =e => віһ20 -l mm , where о 


is defined by sine =R,/R. Eo that we are only considering a case 
where these approximations are very good. 

Let the direction from С to the origin be denoted by u", у", 
м", and let п =u u"+ v у" + ww" be the cosine of the angle between 
this direction and the incident direction. Finally, let E" be that 
new energy which would result from a scattering from the direction 


u, v, м to the direction ц", v", м". 
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We proceed in two essentially different ways according аз this 
energy E" is below or above Ea’ If below, we allow the colliding 
neutron of weight W to scatter isotropically in the laboratory system 
into a new direction u', v', w'. If the new line of flight cuts the 
detector, or if it falls outside of the solid angle subtended at C by 
the detector and corresponds to a scattered energy below cutoff, we de- 
posit the weight W in Lg. Оазе, we follow the scattered neutron 
further to escape or its next collision by the usual (В) routine. It 
will be seen that this is the orthodox way of treating a collision, 
except that we assume that directions within the very small detector 
solid angle have the same behavior as the direction ц", v", м" insofar 
as resultant energy is concerned. 

However, in case the energy Е" corresponding to а", v", м" is 
above cutoff, ме deterministically add a weight Wft to the category Dj» 
which records neutrons impinging on the detector with energies in the 
group containing the energy Е". Here f has its assigned meaning апа 
t - exp {- (в - R,)/a (2) is the transmission for neutrons of energy 
E" in the direction toward the center of the detector. We then allow an 
isotropic scattering into the direction u', v', w'. If this direction 
falls within the counter solid angle, we force a first collision of weight 
Wil - +) at x+ 0"), у + у"), z + w'Àl, where й= - Ал (Е")а [1 - г(1 - t)]. 
(Cf. Chapter III, 85.) If the scattered direction u'v'w' is outside 
this solid angle, we determine the corresponding new energy E' and 
deposit weight W in Ls if E' is below the cutoff E. or follow 


G 


the weight М further in the usual (В) routine if E' is above Ens 
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In а large number М of collisions of the second kind (Е" > Е) 
the expected number hitting the detector will be NWft (with the proper 
energy distribution), as it should be. Moreover (Nf) neutrons will 
scatter into the detector direction, on the average, resulting ina 
total weight of (Nf) W (1 - t) = (NW) f (1 - t) having first collisions 
in this direction, and these are distributed spatially in the correct 
exponential distribution. Finally,an expected number М№(1 - f) of 
neutrons will scatter outside the detector solid angle, with total weight 
N(1 - f)» W = (NW)(1 - f), and these are correctly processed in their 
subsequent history. 

The method has the great virtue of deterministically contributing 
а correct positive weight to the detector on every collision of the second 
kind. (Collisions of the first kind (E"< Ea) cannot do so physically. ) 

A flow diagram covering the method is given in Fig. lla. For the 
и formula, one may refer to §6 of the present chapter. The computational 


parameter Q is set to zero at the source (o). 
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Fig. 44a 
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18. Remarks on thermal neutrons. In all neutron problems with 
which we are concerned, the energy range of neutrons actually followed 
is bounded from zero by some minimum energy E> O. This bound may be 
some stipulated energy E. above the mean thermal energy of the medium, 
neutrons with energies falling below this upon any collision being of no 
interest in the particular problem. The energy E. then serves as the 
lower bound Ес of the lowest energy group for which cross sections аге 
stored. 

However, when neutrons are to be followed down to the mean thermal 
energy of the medium, this energy serves as the minimum Е", and several 
remarks are in order. All cross section considerations and scattering 
formulas have been based on the assumption of target nuclei which are at 
rest in the laboratory system. It is clear that our procedure is un- 
justified in that part of the energy range approaching the thermal energy 
in the case of light elements. Moreover, if the mean energy of the medium 
is in the range below the molecular binding energies involved, simple 
nuclear cross sections may no longer be applicable, (11) 

If neutrons reaching "thermal energy" are dropped, we again have 
Е = Eo - above, neutrons with energies below Е" being thrown into 
a counter Lp for losses to cutoff. 


When the problem necessitates actually following neutrons which 


have reached the thermal energy range, we require а lowermost energy 


(11) 


Theory, D. Van Nostrand Company, Inc., Princeton, N. J., 1952, 
on thermal neutron cross sections. 
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group for such "thermal neutrons," a suitable "average free path," and 
differential scattering laws for this group, аз well as probabilities for 
different types of reactions. It is impractical to deal with the actual 
case of a distribution of neutron energies in the thermal range, and all 
methods assume neutrons within this group have а fixed energy and upon 
elastic collision retain this energy. For heavy elements, isotropy in 
the laboratory system seems to be a valid assumption for elastic scatter- 
ing. However, the choice of proper averages for the other nuclear con- 
stants referred to depends on the elements involved and on complicated 
questions of the actual energy distributions obtaining, which are non- 


Maxwellian in media with strong thermal capture. 


19. Remark on determination of photon sources. Collisions of 
neutrons with particular types of nuclei may result in the production of 
у -rays; radiative capture and certain (п-п) inelastic processes are 
examples. One of the applications of Monte Carlo to neutronics problems 
which is becoming of increasing importance is the determination of 
photon production by neutrons in shields. The latter information may be 
used in turn for the input energy and spatial distribution in the Monte 
Carlo treatment of photon diffusion, discussed in the following chapter. 

Insofar as the neutronics of such shielding media is concerned, one 
need only record as they occur the number of radiative captures, and the 
number of (n-n)radiative collisions caused by incident neutrons of energy 


group g, in each of а set of spatial zones. 
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This, together with a knowledge of the 7 energies produced in 
capture and the production cross sections for 7's by neutrons of energy 
group g in (п-п) reactions, yields the desired distributions for the 


source in the related photon problem. 
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CHAPTER VI 


PHOTON COLLISIONS 


1. Introduction. The procedures discussed in the present 
manual are adapted to problems involving the scattering of photons on 
electrons and nuclei insofar as photons may be considered as "particles" 
subject to mean free path and differential scattering laws. It is for 
this reason that we speak of particles in the text. It is only in the 
collision routine (7) that the physical character of the particles 
need be distinguished. We have discussed many of the collision proc- 
esses for neutrons with nuclei in Chapter V. We now turn to problems 
involving photons. The principal types of photon-electron collisions 
are: (1) Compton scattering, which is the strict analogue of elastic 
Scattering of neutrons on nuclei; (2) pair-production; (3) рһо+о- 
electric effect. These we proceed to discuss from the point of view 


of the Monte Carlo technique. 


2. Basic concepts and constants. We shall use 


с = 2.99776 x 1020 cm sec! for the velocity of light in vacuo 


-1h3- 


Google 


h = 6.624 x 10-27 erg вес for Planck's constant 


28 


ш = 9.10658 x 10 ~~ gm for the rest-mass of the electron, and 


е = 4.8025 x 10710 


esu for the electronic charge 

A photon is characterized (for the purposes of this chapter) by 
its energy є (ergs). To such a photon is ascribed a frequency v (вес?) 
by means of the relation є = Ви and а wave-length А (сш) by the equa- 
tion Лу = с. The "equivalent mass" m (gm) of such a photon is defined 
by me” = €. The speed of all photons in vacuum being с, we may 
assign to a photon а momentum (vector) mV, where V is the velocity 
vector of the photon in the laboratory system. Then 

Ivl = с, and |шу| = mc = hv/c 

An electron is characterized for our purposes by its charge 
e (esu), rest-mass m, (gn), and its velocity vector V. Its mass is 
then m = m/y(1i -B*) where B= |Vl /c. The momentum vector of 
the electron is mV and its total energy me”. 

It is customary to express photon energies by means of a dimen- 
sionless parameter E = hv/n c*, which gives the ratio of photon energy 
to the rest-mass energy of the electron, namely, m c^ = ‚81837 x 10-6 erg 
or .51083 Mev, recalling that 1 Mev is 1.60203 x 10-6 erg (Chapter II,§3). 


This parameter E is the one we adopt for photon energy in the remainder 


of the chapter. 


3. Compton collisions. A Compton collision is by definition 


а collision of a photon with an electron (the latter assumed free and 
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at rest in the laboratory system), with preservation of total momentum 
(vector) and total energy. We may, therefore, proceed in much the same 
way as we did in the case of elastic collision between neutron and 
nucleus (Chapter V, $3), except that we do not introduce а center of 


mass. We agree on the notations and relations of the following table: 


Before collision After collision 
Photon Electron Photon Electron 
а 2 с Xue poe = _ 
Mass m = hv/c m,-m n, = hv'/c n, = m/V1 82 
= м/с 
Velocity V4 Vo = 0 Wi Wo 
Speed MiB: [Уз 20 |ui =с [Ч >| = "2 
Momentum mV. mV, nW nW, 
2 2 А 2 2 
Total energy є = һу = mc mc є' = Би' = nic пс 
We have the two conservation laws: 
шүү Malo = нү oN 
2 2-- 
mc + mc = пус + nc 
or equivalently, 
(hv/c*) V. = (hv! /c*) W + n,W (1) 
l l 22 
hv = hp! + c? (n, - ma) (2) 
2 2 


Introducing the energy parameters Е = hv/n c^ and Е' = hu'/m c^, 


we read 
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EV, = Е' + wA -g? (3) 


вв = 1/J ge сы (4) 


remembering that, contrary to our usual practice, the capitals E, E' 
represent scalars. 


From (№) it follows that 

[1+ (2 - £)]* -13-e (z- z) + (Е-Е!) = g£/(1 -в2) (5) 
while from ( 3) we obtain 

BV - 2 EE' MALA V. BW? = н5 /(1 - g?) 
where Vi is the deflection angle from Và to W- Hence, 


2 - 2 EE' cos ш) = (E - Е')© + 2 ER! (1 - cos V,) 


(E? + E! 
aR asp] (6) 

Eliminating 8*/(1 - BÊ) from the two equations (5 ) and (6), 

\2 


2(Е - E') + (E -Е')2 = (E - E')? + 2 EE'(1 - cos VW.) or 


Е - Е' = ЕЕ! (1 - cos у.) 


Thus we have finally the relation 


E' = E/[1 + E(1 - a)] 
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where we use а = cos у, as usual for the laboratory deflection cosine 
of the incident particle. This is the analogue of the neutronics rela- 
tion E' = E(S + Tu) of Chapter V, §3. 
It is noted that, as in the case of neutrons scattering on 
nuclei heavier than hydrogen, there is a positive minimal energy © 
Е! 
m 


fea Е/1 + 2E for the scattered photon, attained for a = сову = -1. 


As a side remark, observe that the relation 
Е -Е' = EF (1 - cos d ) 


leads &t once to the familiar wave-length relation 
if we first write 


and use the relations E = hv/m сё and v^- с for both primed and 
unprimed variables E,V,^. The "Compton wave length" ть is defined 
to be h/m с = „023263 А (1 Ё = engstrom unit = 10-8 cm), and the 


relation is sometimes written 


у. 
ДА = 2a ва? (51) 
О 2 


We have seen ром a = созу determines E', the resulting energy 
of the photon. We have still to show how the kinetic energy and 
defection cosine cos V 2 of the scattered electron are found in terms 


of E and E'. 
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То be sure, equation (2) answers the first question, since 
c^(n, - m) = hv - hv! = m c^(E - E'); thus the kinetic energy of the 


scattered electron is 
Е. = .51083 (Е - E') Mev 


Ме use equation (3) to determine cos ш. ав follows: 


сов W, = МУ, / LA |у, | = Уз - 6? v (EV) - E'W, ) мос 
Vi -в2 (Ev? zB |у. | КА cos V, ) мос 

V. -в? (к -к'а)/в 

(Е - к'а)/(в / JA. - 5?) 


where &/V1 - в? із known from either equation (5) or (6). 


Ш 


Summarizing the results of this section, we have 


photon deflection cosine а = сов Vi 


scattered photon energy Е' = E/ [1 + E(1 -a)] (units of m c^) 


recoil electron К.Е. E. 


.51083 (E - E') Mev 


electron deflection COS Фо = (Е - E'a)/( в/ - В 2, 


в /}/\ -в2 = Vig - в')2 «2 - £y 


As an example, suppose an 8 Mev photon scatters &t ап angle of 


v, = 60°. Then 
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1 


а 

Е = 15.661 (8 Mev) 

E' = 1.7735 (.90596 Mev) 
Е = 7.09404 Mev 


в/ |р - 8? = 14.854 


4. The Klein-Nishina differential cross section. (12) This is 
the differential cross section e«( 2) in om“/steradian which governs 
the distribution of a = cos v. for the scattered photon in the Compton 
effect, and is the photon analogue of the neutron cross section of 
Chapter V, §§4 and 6. It is defined by the formula 


2 2 2 
TONO 2—22 { . . Ez -af lio (7) 
Е Ё 4) у (1 + а2) [1 + Е(1 -а)] 


where г. = e^ /n c? = .28183 x 10712 om. 


Now we might proceed just аз in the case of elastic collision 
to determine а by (say) von Neumann's device (Chapter 1,85) апа then 
E' from E and а using the relation 
(12 or ап account of this and related functions see В. Latter and 


Н. Kahn, Gamma-Ray Absorption Coefficients, Project RAND, R-170 (1949) 
and the National Bureau of Standards Circular 542, Graphs of the Compton 


Energy-Angle Relationship and the Klein-Nishina Formula from 10 Kev to 
500 Mev. 
(13) mis is the classical Thomson radius of the, electron which enters 
into the electron cross-section formula (8/3)nr* = .665 x 10° cm? 
(= .6654 barns). 
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Е' = Е/[1+Е(1 - a)] (8) 


of the preceding section. In view of the complexity of (7) it is 
fortunate that a simpler alternative exists. We describe this in the 


following section. 


5. The photon energy distribution and Compton cross section. 
In order to prepare for the alternative method, we first determine the 


differential cross section б(Е') from the relation 
О ! "= ==» e 
Q(E')aE' = o,( 9)49 e (9) * 2 т(да) 


Using the preceding equation (8) in the form 


(9) 


we obtain the result 


2 
nr 2 
7 (mt) = © E' 2, i 2. 1 1 
ate" в2 (Ё "(+ 3)-Ё Е 2)(3.) +} 


We may obtain easily from this the total cross section for 


Compton scattering as follows: 


Е pm 
e „(Compton ) = ] c ,(E' )dE' 
E 


1 + 2E 
2 1+Е 2 E^ - of - 2 
uer. uode pr cni Ов) 
(1+2) Е 2Е 


Google 


which is graphed in the National Bureau of Standards Circular 52 

( 1с ) and may be used to compute free paths ae = 1/29. (Compton) 
for & suitable set of energy groups (assuming No is the numerical 
density of free electrons and no other processes exist). 


We will actually use the Monte Carlo procedure 
E m 
г = P,(E") = о (E! )аЕ! e, (Compton) 


to determine E' from E and random number г, since a sufficiently 


(14) 


accurate fit for the inverse function is given by 


; E 
E а 
1+ 4х + (8E - г 


where # = Е/(1 + .5625 E), and E $a (—2 Mev). Addition of & term 


Е(Е-А) г (1 - x)? 


yields a reasonably good fit on the range }4< E < 10. 

The equation (5) above then permits determination of a from 
Е and E'. 

Thus a Compton collision with incident energy E not exceeding 


5 Mev may be satisfactorily handled by the flow diagram of Fig. 45. 


ee eee oe eee НОЖИ 


ny Ee Oe Ee Ee Ee 


and Energies, LAMS-1199, Los Alamos Scientific Laboratory, 1953, for 
goodness of fit. 
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Fig. 45 
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6. Photoelectric effect and pair production. We have discussed 
the Compton collision for photons with electrons and noted that it 
corresponds to the elastic collision in neutronics. If a medium con- 
sists of elements A, B,..., the contribution of the Compton effect 
to the "total cross section" 2 (cf. Chapter 111,82) is Ne, (Compton), 


where N - ААА + Nan +... is the total numerical electron density. 
3 


Here М denotes the number of atoms of element А рег сш’ and Za 


is the atomic number (= number of electrons per atom) of element A. 
In determining the "total cross section" all collision processes 


must be taken into account, and one must consider in this connection the 


further contributions мо (ре) +... 


effect and мот, (pp) +... =} (pp) of pair production. Cross sections 


E 
for these processes for elements with 4< Z < 92 and .10 < E Á 20 (m cî) 


У (pe) of the photoelectric 
E 


may be found in the RAND report R-170 referred to before. Thus the 


free path determining collision positions is given by 


A = ME 


where 


Х = Na (Compton) +5 (ре) +}; (рр) 
E E 


In the present report we regard pair production and photoelectric 
effect as absorptions, and have therefore, upon any collision, an 
absorption probability [Е - No.(Compton)]/Z , which may be dealt 


with by the alternative methods discussed in Chapter II, §2. 
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It must be remembered that secondary photons (X-rays) may result 
from the photoelectric effect in heavy elements, and that the electron 
and positron born in pair акыйнек csoduce bremsstrahlung while being 
Slowed down by ionization of the medium. Moreover, the inverse photo- 
electric effect may come into play, and the positron is finally annihi- 
lated with an electron to produce still further y-rays. If the energy 
cutoff of the problem is not sufficiently high to exclude consideration 
of these secondary photons, one must include in the output of the orig- 
inal problem sufficient information to determine source distributions 
for such secondary radiation. We do not deal with this case in the 


present manual. 
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CHAPTER VII 


DIRECTION PARAMETERS AFTER COLLISION 


1. Introduction. The purpose of the present chapter is to 
develop formulas for the final direction parameters of a particle after 
scattering through an angle MI of cosine a, in the laboratory system, 
from an incident direction u, у, м (ог м alone). While these formu- 
las are somewhat long, they may be derived from the simple principle of 
elementary complex variables which states that (x + iy)(cos O + i sin Ө) 


is & complex number whose vector is rotated through an angle O from that 


of x + iy. 


2. Formulas for the final direction cosines. Let u = cos а, v = 
cos B, w = cos у be the direction cosines of the incident line of flight. 
Consider (u, v, м) ав а point on the unit sphere ee + y? + w^ = 1 in 


direction space U, V, М (cf. Fig. 46). Letting Р, $ be the polar 


coordinates in the U, V plane of the point (u, v, о) we see that 
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р = sin y = yj. 
cos $ 


iu 
£l 
NS 
hd 


sin 9 у/ р 


Our first objective is to derive formulas for а particle rotation 
of U,V,W space into itself which takes the point (0,0,1) into (u,v,w); 
namely, that resulting from iteration of two simple rotations: the 
first about the V-axis through the angle у, followed by а second 


rotation through the angle $ about the W-axis. 


Fig. 46 
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We have for the first rotation 


w' + iu' = (w + iu)(cos У + i sin ») 
} (1) 
v'= у 
and for the second 
u" + iv" = (u' + iv')(cos Q + i sin 9) 
| (2) 
у" = м! 
Separating real and imaginary parts, we read 
u' = u cos У + w Sin 7 
pus v | (3) 
w' = -u sin У + W cos У 
ч" = u' cos $ - у! sin ў 
v" = u' sin ў + v' cos Г (4) 
м" = LA 
Substitution of (3) into (h) yields 
ц" = u соз У сов $ - v sin +w sin? сов Q- 
у" = u cos y sin + v cos @ + w sin” sing $ (5) 
м" = y 


-u sin * + м COS 


Replacing the functions of 7, ў ъу their values in terms of 


TE va _ e ‚ ч, у, and м, we obtain from (5) 
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с | 


ц" (имо - уу)/ р + м 


(6) 


< | 


wv¥tvu)/e + м 


< 
i 

=` 
£ 


Under this rotation (6) it is evident that (u,v,w) = (O, O, 1) goes over 
into (u", v", м") = (u, у, м). 
Now we know that the direction (и, V, W) of the deflected line of 


flight makes an angle W., of cosine a with the direction (u, v, w) 


l 
of the incident line of flight. It is clear that this restricts the 
former to a cone of directions of opening va about the latter direction, 
and the directions on this cone are all equally likely. We adopt the 
following arbitrary convention for fixing the deflected direction; 

namely, we always consider first & cone of opening Wy about the 

cos 6, sin w 


Weaxis OW, and а point Р = (sin y sin 6, cos v.) 


1 L 


located on this cone and the unit sphere, determined by an azimuthal 


angle of 6 uniformly distributed on -xm $ ё 5 т (Fig. 47). 
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Fig. 47 


We then agree that the final direction (u, v, w) is the image 
under rotation (6) of the point P so constructed. It should be 
clear that this is an appropriate convention for our purpose. 


For simplicity we adopt the further notation 
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а = соз vai 


MN" _ 2 < < 
b = siny = (1-а) OS wm 
с = cos 6 
а = 5118 = (seni )V(1 - ét) -п<6 <т 


In this notation, the point Р = (bc, bd, а) = (и, у, м) goes over into 


ч = (bewu-bdv)/p + au 
У = (осму + bdu)/ р + ау 
w = -bcp + aw 


Summarizing the results of the present section, we have 


(bewu - vav)/ Aa = w^) + -— 


v! = (осму + bdu)/ ya - w^) + av (7) 
wt = -be e Е w^) + ам 


where 


u, у, м аге direction cosines of the incident line of flight 
ut, v!, w! are direction cosines of the deflected line of flight (lab.) 
а = cos W,,where Фу is the angle between the two 


b =? = а” 


с = сов ё, where -rm $a € л uniformly 


а = (sgn ё ) Ёл” 
Note that these formulas should not be used if | | is too close to 


unity. They are then poor computationally and indeterminate at MB Ia 


In such a case it is better to by-pass the rotation (6) and use for 
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ut, у*, м! the coordinates 


ut = He 
v! = bà 
у! = ач 


where а, 5, с, а have the definitions above. 


3. Subroutine for the final direction cosines. We incorporate 
these results into the (ô) routine of Fig. 48. The power 27] used to 
determine whether the incident direction should be considered vertical 
may well depend on the accuracy required and on the number of "bits" 


afforded by the particular machine being used. 


-161- 


Google 


Fig. 48 
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hl. Final direction м in slab or spherically symmetric case. 
The present section is concerned with determination of the final direction 
cosine м! = cos у! of the scattered line of flight as it exists at the 
point of collision, in terms of the cosine a cf the laboratory angle v. 
determined in (У), аз required in both slab geometry а Brlerfcsiis 
symmetric geometry. (Cf. Chapter II, 82, and Endet es ТУ, $8.) We recall 
that, in the spherical case, the cosine м = cos y of the angle у which the 
incident line of flight makes with the radius vector as it exists at the 
point of collision has already been prepared (Chapter IV, $3) and stored as 


w before entry at ($). 


It is immediately evident that the м! formula of equations (7), in 


у! = e Vi "mem 


applies to the slab case, where 8, b, and c have the definitions given there. 


$2 of the present chapter, 


Note that the formula is independent of the u, v used in the derivation, as 
it should be. 

But it will also be clear that the identical formula applies to the 
spherical case, if we imagine for the moment that OW represents the outwardly 
directed radius vector direction and u, v, w the direction of the incident 
line of flight relative to the radius vector, the u,v being immaterial. 

In both applications, it will be observed that the angle 6 involved 
in c -cosó may be limited to the range O < g < л, since cos -6 = cos ў = с 


is the only function of 4 appearing in the м! formula. 
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Fig. 49 


In both cases, therefore, we have the simple routine of Fig. 19. 


5. Scattering isotropic in the laboratory system. Іп any case 
where particles emerge from a collision in a distribution isotropic in 
the laboratory system, it is patently foolish to determine a = er - 1 
for the deflection cosine relative to the line of flight and then determine 
the final direction cosines by means of the preceding two sections. It 
is simpler to Heri es direction of the deflected line of flight by the 
same method we used to set up an isotropic source in Chapter II, $58, c, 
ignoring the incident direction completely. Thus the ( 8) ЕРЕ ЕИ referred 
to is simply that of Fig. 10, with exit to ( B} for the u, у, м case,and, 


still more simply, boxes г and м = 2r - 1 for the slab and spherical case. 
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CHAPTER УТТТ 


TERMINAL CLASSIFICATION 


1. Introduction. We have already indicated various terminal 
events in the history of a particle, e.g., transmission, capture, loss 
to energy cutoff, and so on. We must still consider some of the many 
ways in which it may be desired to classify a particle which escapes 
the system. This is the general function of the (€) routines pre- 
viously referred to. Such routines invariably exit to (а) for the 
introduction of а new source particle. Of these many kinds of classi- 
fication, we can hope to give only a brief indication, since the demands 
of physicists on this score are frequently involved and exacting. Indeed, 
it is the ability of Monte Carlo methods to provide answers to the most 
intricate questions of this kind which makes them an indispensable tool 


in design. 


2. Classification of escapes on number of collisions. As the 


simplest example, we may refer to a problem whose purpose it is to 
determine the distribution of escaping neutrons with respect to the 
number of collisions suffered within the system. Suppose we agree to 


keep storage registers Ny» v = l, 2, ..., 10 for the total weight of 
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neutrons escaping with и < 10 collisions, and м] Рог those escaping 
with v 2 11. We suppose that the forced first collision device has been 
used so that the transmission counter Т records the weight escaping 
with у= О automatically (cf. Ch. III, 84). We may, therefore, follow 


the routine of Fig. 50. 


© 
о-о 
© 


Fig. 50 
3. Energy and angle distributions of escape. Consider a problem 
with a source direction u = 0, у = 0, м = 1 in which one desires the 
correlated distribution of escapes in energy and angle with the source 
direction. There must be provided in permanent storage a suitable set 


of cosines C, » C, >... > С 


1 2 I -l and energy bounds Е >>... 5€ 


Н 
= Ес, which need not be the same as those used for other purposes, while 
in dynamic storage, we reserve JH positions Ns oh for the total weights 
escaping in the corresponding categories. Іп such a problem, an escaping 
neutron enters (є) with a known weight W, direction cosine w, and 


energy E, and is classified as indicated in Fig. 51. 
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Fig. 51 


Such a classification in angle is appropriate in a case like that 
of Fig. 52, where a detector band is at essentially infinite distance 
from a scattering medium symmetric about the Z-axis. Here infinite 
distance means that all particles escaping from any point of the 
medium surface in the same direction hit the band in the same angular 
zone. Then the number Nan referred to is the number of h-group 
escapes in the solid angle 2т (С, - Cs) so that N, p/e" (С, - С 1) 
is the number in this category escaping рег steradian, and the number 
hitting a detector of given area on the circle indicated in the figure 
can be predicted. 

If the source direction is again u = 0, у = 0, w=1, and the 
scattering medium is symmetric about the Z-axis, ensuring symmetry of 


escape about this axis, but the detector is not at infinity in the 
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uniform 


detector band 


Fig. 52 
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sense referred to, it is necessary to take into account the distance 
ав from center О to the detector band as indicated in Fig. 53 ара 
Fig. 54. In the latter, the exit leads to the previous classification 
routine of Fig. 51 and ессе to (а). At entry to the routine of 
Fig. 54, the x,y,z parameters are those of the last point of de- 
parture before escape. Again the results may be normed to number per 


steradian and interpreted with reference to the detector band indicated. 


detector band 


Fig. 53 


Google 


То Fig. 51 


Fig. 54 


The preceding discussion is limited to systems symmetrically 
distributed about the source (+7) direction. Problems of this kind 
have the great advantage, from the Monte Carlo standpoint, that no 
escaping particles are lost to classification. 

Unfortunately, there are experimental considerations which, in 
many cases, dictate a geometry in which such symmetry is lacking. 
Consider, for example, the situation in Fig. 55, in which a source 
with direction u = О, у = 0, w= 1 impinges on the lateral surface 
of a cylinder with axis on the X-axis. It is desired to count the 
numbers "s hitting a coaxial band zone j upon escape, the de- 


tector band having radius аз and height h. Since the cylinder 
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non-uniform 


— — => — J — —- — — — — 


180° 


Fig. 55 
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lacks symmetry about the Z-axis, it is no longer possible to classify 
escaping particles according to the cosine м of the angle у of the 
line from origin to the point of intersection of the line of flight 
with the sphere of radius da» since the escapes with this direction 
are non-uniform about the Z-axis. Without rather complicated devices 
for prejudicing scattering in the band direction, which we shall not 
discuss, all one can do is to submit to running very much larger 
samples in order that the rather small solid angle actually subtended 
by the band shall receive а sufficient number of escapes to render the 
М, statistics reliable. This is а real difficulty, and it would be 
of great value for experimentalists contemplating the use of Monte Carlo 
methods for thick target corrections to consider the feasibility of the 
symmetric geometry in the experimental setup. 

A suitable routine for the case cited, without any special 
importance sampling devices, is given in Fig. 56, the exit again re- 
ferring to the routine of Fig. 51. The category La refers to all 


escapes failing to hit the band. In computing the distance t to the 


(infinite) cylindrical surface y^ + a = one avoids the 1 - Ts = 


catastrophe, аз indicated in Ch. IV, 86 (end of section). 


О 


Consider finally а case where it is desired to classify а 


particle escaping with direction u,v,w from the lateral surface of 


a cylinder у + та = Ro according to the angle measured from the 


normal to the surface. If (X,Y, 92 is the point of escape, deter- 


a) 


mined in the usual manner from u,v,w,R, and the position x,y,z of 
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To Fig. 51 


Fig. 56 
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last departure, then the direction of the normal at this point is given 


by 


and the cosine w of the angle of escape with the normal is 


X,u + y,v +2, • О 
t t t 
н uud АА 
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СНАРТЕВ ТХ 
REMARKS ОМ COMPUTATION 


1. Scaling. If the machine employed lacks the "floating decimal" 
feature or if considerations of machine speed exclude its use, all 
computations indicated in the flow diagram at each step must be planned 
in advance to remain on the interval -1< x <1. The formulas and flow 
diagrams of the present text are unscaled and are expressed in the usual 
physical units. The scaling procedure in Monte Carlo problems usually 
causes no difficulties, and may be left to the imagination, with the 
reminder that one must be careful to avoid the loss of accuracy which 
attends over-scaling. Remarks are made about scaling at various points 
in the present chapter as it enters into special subroutines, such as 
exp( - x/y) where x may exceed y, in the Qn x routine where x 


2-23 


may be аз small as » and in the cosine routine for - т <x < л, 


2. Debugging. When a problem is ready for actual running by 
the machine, some method of detecting the inevitable human errors which 


may enter at all phases of the coding must be employed to assure all 
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concerned that the machine will perform exactly the routine intended. 
It goes without saying that all permanent storage constants should be 
printed out at various intervals and checked against the original list, 
certainly at the beginning and end of the problem. Moreover, it is 
sometimes possible to build in a sum check of all code orders and 
permanent constants which will automatically detect any changes due to 
electronic misbehavior once the machine is running. We refer here 
rather to the question of ensuring that the code itself is correct. 
This question presents rather greater difficulties in Monte Carlo 
problems than in other types of computation because of the complexity 
of. decisions involved. 

We have found that the most convincing guarantee of faithfulness 
of code to flow diagram consists in the preparation by hand of a 


deterministic sequence of numbers т), P5, T designed to process 


3c 
automatically & hand-picked set of particles which are so chosen that 
every logically possible path through the flow diagram is traversed at 
least once with non-trivial parameters. A special debug routine may be 
used for this purpose which instructs the machine to follow the code 
precisely except that, upon call for the next random number, the 
random number routine is by-passed, and the prepared list is consulted 
for the next number r of the hand-picked sequence, At every list 
consultation the machine prints out all significant variables of the 


problem, the printout then being compared with the hand-computed 


particle histories. 
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This is a laborious process іп the hand=computing phase, but is 
well worth doing, even though, in complicated problems, the processing 
of 50 particles ог so may be necessary. 

3. Special subroutines. We include a number of special sub- 
routines for some of the functions commonly occurring in Monte Carlo 
calculations. No claim is made that these are the best available; they 
are only given for completeness in case the reader knows of no better 
one. 

a. A random number routine. We have used in all problems the 
sequence гу, ro, ETE of 38 bit diadic numbers, generated by the 


algorithm defining Tr to be that 38 bit number corresponding to the 


+1. 


middle 38 bits of the square of г, where г. is the 38 bit number 


1 
defined by 10 BBB FALDE in the system with base 16, with the decimal 
after the first bit. Thus гу = 0.001, 0000, 1011, 1011, 1011, 1111, 
1010, 0100, 1101, 1118, the final zero being ignored. The sequence 
automatically terminates in zero at about n = 750,000 and has been 
thoroughly tested, not only for the usual statistical features but by 

its actual use in many problems. The Monte Carlo use of the sequence is 
of a very peculiar and unpredictable kind. If one fixes attention on 
any particular random number box of the flow diagram, it appears that the 
numbers actually selected from the sequence for use in this box are а 
subsequence selected from the main sequence by completely inscrutable 


rules depending on the sequence itself. Whenever a probability of some 


terminal event (e.g., transmission) was known, the corresponding N, /N 
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for that category agreed even better than one might expect it to on the 
basis of a truly random sequence. 

Ъ. Shifted random numbers. Ап average problem may involve 
20,000 source particles, each with perhaps five collisions, and each 
collision may involvethree or four random numbers. It is clear that 
one may expect to exhaust the random number sequence in many problems 
if it is used in this simple way. it is true that the same sequence 
may safely be used over again without repetition of results provided 
the initial ri is not used &t the identical place it was first called 
in the flow diagram, and indeed, such re-use of the sequence may be 
resorted to. 

We may mention however that the usual problems do not require 
anything like 38 bit resolution, and the number of random numbers 
available may be increased by a factor 2 or even № by а shifting 
routine, which is incidentally faster than squaring and so helps to 
reduce machine time. This consists in using in sequence not only T? 

E a2, ©] 
but also the fractional parts of 2 г, 2 (2 т), etc.,as random 


numbers, before squaring т, to obtain the next г For example, 


n+l’ 
if three random numbers are to be obtained from r,» one might use the 
fractional parts of ra? a'r and г (г бк). Ihe sequences ob- 


tained by such extensions of the basic routine have also been tested 


and proved satisfactory. 
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с. А logarithm routine. From the usual series 


fa +х) = х - x^ /2 + х? /3 = xt fy +... | x | <1 
one obtains 
fa ox) шеи еа Bosch does Ир Ix] «a 
and hence 
fo (==) = - 2 {x + к3/з +525 + sss) |х | < 1 
C Lex 


which, under the transformation у = л becomes 


E _1\3 742 
ЕССЕ 


Convergence of the latter series is rapid enough for the range 
= < у < 1 to permit use of a moderate number of terms; five should 
suffice for most purposes. 

Now suppose t is a number on the range 2-23 < t < 1 for 
which fh t is to be computed. We determine that integer n for 


which 


2-23 < g-(m1) < в , $n 


IA 
- 
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Then ie 27 & =n <1 and j^ ё = - n fh 2+ 41, where 
fh п can be obtained from the series. 

Now - 4, а ue 23 Ёл 2 < 23(.694) < 16, so that |2-* Rn t | < 1 
for & on the range 2723 < Ёё < 1. Hence а routine suitable for а 
scaled problem is given by Fig. 57. 

Our only use for the Rn function has been in collision- 


distance formulas such as 
Ё - "T . 2" х (27 фт) 


Limitation of random numbers entering this formula to the range г 2, 2723 
means that we force a collision within 16 free paths. This is inaccurate 
to the extent that е © is the chance of collision at a still further 
distance, which is an acceptable error in most cases. Іп problems in- 
volving systems so large that such errors are significant, probably some 
type of importance sampling will be used in place of the simple formula 
above. In any case it is clear how to modify the routine for higher 


powers, 2723 


being quite arbitrary. 
For the routine as given, it is clear that if all distances (сш) 
are scaled down by a factor exceeding 16 times the greatest free path 


A occurring in the problem, one is safe so far as this particular 


formula is concerned. 
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а. The exponential ехр(- х/у). It is frequently necessary to 
compute the value of exp(- x/y), where O< x<1, О<у< 1. Even 
though x and y are scaled, the quotient x/y may exceed unity and 
cannot be dealt with directly by fixed decimal point machines. There 


(15) 


are excellent polynomial approximations to exp(-p) for O<p 5 £n 2, 


so that ме may define the integer п by 
n La 2 < хју < (1+1) fh 2 


and р оп the prescribed range by Р = (x/y) - n Ln 2. We shall 


therefore have 


expl- x/y) = 277 е7” 


-p А | 
> а. 
where е а, +p ( 1 + ра.) with 


a = l 


- 9664279798 


w 
li 


@ 
| 


> = 3535763631 


is usually sufficiently accurate for our purposes. 


[GM ENTER PEN eee Е E 
15/8. Carlson and M. Goldstein ; Rational Approximation of Functions," 


Los Alamos Scientific Laboratory, LA-1943, 1955. Gives polynomials 
of various degrees with bounds on the errors involved. 
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We may therefore proceed аз in Fig. 58. 

е. А cosine routine. We indicate a scaled routine for computing 
2 cos 6, where 6 is to be chosen uniformly on -7 < ё < 7, as 
required by Ch. VII, 8 3. It is based on the following series of 


formulas 


Input 0zx«1,0«y«1 


Stored constants a ,a_,a „2n 2 
о 1 2 


"Кш + + 
а, pla, ра.) 


-р 


ехр-х/у) = 2 Пе 


Fig. 58 
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the final formula arising from 


cos ё = cos Ш ё! UT E б' - 1 


2 
-l + 2 Ё - 2 ein e| 


3885488 


It may perhaps be mentioned here that for scaled machines the 
following procedure is somewhat shorter than that obtained by applying 
the preceding method to Fig. 19, the routine for slab and spherical 
final direction w'. In place of choosing 6 uniformly on OS 6 < m, 
we may introduce the angle 8 = 6 - л [2 ‚ choosing the latter uniformly 


od 
on -7/2 < 5 < 1/2. Thus we may proceed as follows: 
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S = sin g' as above 

Ad 
с = cos б = sing = sin 26' = 2 sin 6' cos ё! 
w! = ам - bc l-w 


= aw-2slf(1- a j(i = s^) - м”) 


We should therefore have in place of Fig. 49 the routine of 


Fig. 59. In this way we avoid the square root for b = yi = а? апа 


make use of only one double angle formula instead of two. 


Google 


ц. А Monte Carlo device for yr . There are many ingenious 


devices for avoiding the computation of the value of various functions 


(16) 


with arguments involving a random number. As a simple example, the 


process 


involved in throwing for the radius of a neutron in a source distributed 
uniformly in area on a disk of radius Ri may be replaced by the fol- 


lowing: 


E г т R = В шах(ё, n) 


That is to say, опе may use in place of the square root of one 
random number the larger of any two successive ones, for the equation 


é = yr results from the standard Monte Carlo relation 


x 
г = | СЁ de 
о 


and is thus equivalent to choosing & оп the interval (x, х + dx) with 
frequency 2x dx. But the alternative of throwing points (2,95) into 
the unit square О < &<1, О < п < 1 uniformly in area and choosing 


the larger of the coordinates automatically determines an x = max(t,m) 


(16 ), 


e Je М. Butler's paper in Symposium on Monte Carlo Methods, 
John Wiley & Sons, Inc., New York, 1956. 
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on (Ё,Ё+ dé) with this frequency, as is apparent from Fig. 59a. 


5. A Monte Carlo device for the cosine of ап equidistributed 
angie‘ 17 ) In various places we have needed the cosine and sine of ап 
angle ф  equidistributed on the interval -mr < ф & я, and have hereto- 
fore given a rather cumbersome but straightforward method involving 
series, square root, and, in case of fixed decimal requirements, a 
multiple angle formula. We give in Fig 59Ъ an alternative method, 
easily scaled if necessary and avoiding series and square root routines. 

One observes that the desired determination of с = соз ф and 
d = sin$ 15 equivalent to choosing а point (c,d) uniformly on the unit 
circle. This in turn is tantamount to throwing points (#,7) uniformly 
in the square -1< x, y $1, rejecting those lying outside the unit 
circle хе + 26 = l, End taking for c and d the values $5 Ый n? 
and n Wes пе, for the retained points (6, л) are uniformly distri- 
buted in area in the unit circle and hence their projections (c,d) on the 
unit circle are uniformly distributed also. But these square roots may 
be avoided by limiting the selection procedure to the first quadrant and 
using double angle formulas to obtain (c,d) uniformly on the upper half- 
circle. Finally an additional random number can be used to change the 
sign of d with probability 1/2. The efficiency of the retention of 
random number pairs is т/і. 


Commerce, National Bureau of Standards, Applied Mathematics Series #12, 
Monte Carlo Method, Washington, D. C.e 1951. Cf. footnote (1). 
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22 


Shaded area = 2x dx 


Fig. 59a 


|- HH + H 


Fig. 59b 
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CHAPTER X 


STATISTICAL CONSIDERATIONS 


1. The limit theorem in the Bernoulli case. Suppose that a 
certain experiment can result in c ways Of which a are favorable 


to an event E, while b=c- a result in the event Е, (по% Еу). 


Consider the set of all sequences of М trials of this experiment. 
In this set of en sequences, the number of sequences resulting in 


exactly M occurrences of E is clearly 


l 


where e is the number of combinations of М things taken M ata 


time. Hence, the probability of exactly M occurrences of Ei in a 


sequence of N trials is 


N N M ,N-M,jN N M М-М 
Ру = CM а b /C = CM ра 


И 


where р = а/с and а = b/c are the probabilities of E, and Е, 


respectively. 
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It is meaningful, therefore, to speak of the probability that, 
in a sequence of N trials of this experiment, the number M of 


"successes" E, shall lie between two given bounds; indeed, 


l 
P(a« M« Bg) = РА a 


^. а<М<В Рм 


18 
There is a fundamental probability thecren! which states 
that 
t 2 
Е - > < є) = \ 2 | -u fe du + В 
9 
where t = є Y— and 
2 Ира 
Е 22+ .25 Ip - gl -6/2) Ура 
[ у + а + е 
emNpq si 
2. Application to the terminal ratios. The output of a non- 
multiplicative! 29) Monte Carlo calculation without weights gives the 


До 
16), 


. V. Uspensky, Introduction to Mathematical Probability, McGraw- 
Hill Book Company, Inc., New York, 1937, p. 130. 


(19) 


I.e., not involving fission,(n-2n) reactions, etc., in which а 
particle gives rise to more than one particle. 
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number M out of N source particles which terminate in each of а set 
of all-inclusive, mutually disjoint categories C. If we fix attention 
upon some particular category C, we may regard the processing of a 
source particle as an experiment which has as its outcome either the 
event Ei of termination in this category or E of non-termination 
therein. 

Now in any actual problem, there is & definite upper bound l on 
the number of random numbers needed to process a particle due to the 
existence of cutoffs based on energy, time, weight, number of collisions, 
etc. We may then consider the class of all ? sequences of | random 
numbers, All these sequences are equally likely, and а certain а of 
them determine а history terminating in category C. We may, therefore, 
say that the probability of termination іп C is pz=a/y. 

To be sure, this probability p of E; is unknown; indeed, its 
determination is precisely the object of the problem. We may, however, 
gain some idea of the statistical reliability of the Monte Carlo result 
by tentatively taking for р the value of M/N at — late stage of 
the problem, when the latter ratio appears to have settled down, and 
define а = 1 - (M/N). Then the preceding theory states that the ratio 
of the number of sequences of N trials resulting in a ratio of B 


N 
satisfying the inequality 


Ё > P € 
N P 
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to the totality of all possible sequences of N trials is approximately 
f(t) = erf(t/y2) 


where + = є e and erf(x) is the well-tabulated "error function" 


2 [* у 
erf(x) uri e dx 


for which we include a brief table. 


Fig. 60 
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Ву assigning to є the maximum error to be tolerated and to N 
the number of source particles processed, one finds from the value of the 
integral the approximate chance of ап error not exceeding є, which should 
be close enough to unity for comfort and which approaches unity with in- 
creasing N. 

As an example, suppose a Monte Carlo run of 50,000 neutrons shows a 
capture of 5000 neutrons in a certain zone, and the probability p of cap- 
ture in this zone is desired with 5% accuracy. Let р = „1, е= .05p = .005. 
Then t = 3.72678 and f(t) = erf(t/y2) = erf(2.635) > .9998. The 
error H in using f(t) for P(lM/N - p| < є) does not exceed .0001. 

This is, of course, far higher probability of safety than one can 
usually hope for. Consider for contrast the extreme case of а capture p 
of about .Ol with & maximum error of 14. Then є = 10^", and for N - 50,000, 
t = .225,and f(t) = erf(t/y 2) = erf(.175) ® .18, with |R| < .018. Here 
it appears that & simple-minded Monte Carlo is quite ineffective. To be 
sure, these requirements аге far more stringent than are usually encoun- 
tered, but such problems do arise, notably in counter design, and one sees 
clearly the necessity for very large samples combined with ingenious de- 


vices for improving statistics by use of weights in such cases. 


3. The central limit theorem. The preceding remarks apply 
strictly only to Monte Carlo procedures of the most straightforward 
kind where no weights аге employed, and а single, non-multiplying source 


particle is followed until it ends its history in some terminal category. 
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It may then be considered аз a physical particle undergoing ап experiment 
insofar as the random numbers employed may be regarded as truly random. 

However, we have mentioned many devices employing weighted par- 
ticles. In such cases, it is clear that the total weight M terminating 
in a given category C after such procedures no longer represents the 
number of successes out of N trials, that is, out of the processing of 
N source particles, each initially of weight 1, and we can no longer 
apply the theorem of 81. There is, however, a well-known generalization 
called the central limit theorem, of which we shall give here a very 
special вава, О най does apply to the weight method. 

Instead of the simple experiment which results in 1 with prob- 
ability p andin O with probability q, let us consider an experiment 
which can result in K different ways, to which we assign definite real 


numbers W Saas Wes with probabilities Ру, —— Рк, respectively, 


1? 


р] + coe + РК being unity. In a single trial of the experiment, therefore, 
we may say that Whe has probability Py» К = l, ..., Ke In such а case, 
we define the mean а = Pp. W and the dispersion b = E pW, = аў” 


= Z р, ie - aÔ, Suppose now that N trials are made of this experiment, 


and М is the sum of N weights so determined. Then our theorem states 


that the probability 


Е - а < є) = erf(t/y2) + р, 
(20) 


J. V. Uspensky ( 1c), page 294; cf. footnote (18). 
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where + =e уг апа Pu О as N > о, The estimation of Р 15 
given in terms of the third moment of the W-distribution. 

Notice that the Bernoulli case is contained in this, since for 
that case, we have р, = Ps Wi = 1: Po = Qs М> = О, а= р. 1+а • 0 = р, 
Ъ = (р. is +q» 0°) ar р - uis = p(l - р) = ра, and the sum M of 
the weights 1 and О recorded in the М trials is simply the number 


of 1's (successes) in this sequence. 


4, Application to problems using weights. Consider now a per- 


fectly general Monte Carlo problem 21 


in which weighted particles may 
be used, each particle leaving the source with weight 1. For simplicity, 
we consider the weights limited to a finite set of numbers, which is, in 
fact, the case in machine computation. Each possible history assigns a 
definite weight W to the particular terminal category C. Moreover, 


it is clear that each such weight М, has a definite probability Py 


k 
of realization upon the trial of а random source particle, namely, the 
probability of & history which assigns W to C. Thus the mean 

а = $ Py М is the expected weight terminating in С рег source раг- 
ticle. If we process N particles in a Monte Carlo problem, tallying 
the weight wit) contributed to category C by each trial particle 
T=1, 2, ..., М, the final sum zw) is the М of the theorem, and 


M/N may be used as an approximation to a. The difficulty in applying 


(21) 


Multiplicative processes are not excluded. When they аге involved, 
we speak of the "expectation" a of termination in category C 
rather than the probability whether weights different from unity 
are used or not. 


-195- 


Google 


the theorem lies in the fact that ме do not ordinarily have at hand ап 
estimate for the dispersion Ъ. In the simple Bernoulli case, we saw 
that b was expressible in terms of p, but, in general, a does not 


k 
as a guess for bo if we had provided in the problem for tallying the 


we 2 Е т\2 2 
determine b. Since b = Хр, М -а , we could use J (W")"/N - (M/N) 
T zl 
т\2 т : : 
(И) аз well аз the (W) for each trial particle T= l, ..., М. 

Note that in case of capture, one squares the total of capture 
contributions from all collisions of the particle; it is not correct to 
simply cumulate the squares of each weight captured. 

In a problem which does not involve multiplication (fission, (n-2n) 
reactions, etc.), the weight procedure is better than that not employing 
weights, at a given stage М, and for a given error €, to the extent that 
the upper limit of integration ү exceeds € vx. the probability in 
the left member being identical, that is, the gain is based on the degree 

А { we 2 2 
to which b is less than ра. Now bz ХР). ха, ра=р-р, апа 
the expectation of category С in such a problem is its probability, 


1.е., а = p. Hence, finally the improvement rests on the fact that 


which is clear since W, < 1 Рог all К. However, the degree of 
improvement can be calculated only in the most trivial examples. It 
is at least clear why the improvement is so great in cases of small 


quantities, since, then, the We are quite small and м. << We 
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It may also be remarked that the theorem of section (3) applies 
perfectly well in finding the expected number a (per source particle) 
of particles terminating in a particular category C, while the theorem 
of $1 does not. The unfortunate feature involved in the application 
of the former theorem is that one has no simple way of knowing the 
dispersion b. Problems of the complexity in practice usually forbid 
complete tabulations of Wig for purely statistical studies. Actually 
one provides the best weight tricks one can think of and tries to judge 
the reliability of an output ratio by its convergence and stability as 


N increases. 


О. [Illustrative examples. We include some extremely trivial 
examples which nevertheless may help to fix the ideas involved in the 
preceding sections. 

(1) Consider the non-multiplicative problem defined by the 


following flow diagram: 


ec: 
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Ме may limit random numbers to be either О or 1. The upper 
bound | on lengths of sequences of random numbers required is 2, and 
the total number У of such sequences 15 22 = №. In fact, we have the 


following results: 


OO implies termination in A 
Ol implies termination in A 
10 implies termination in B 


ll implies termination in C 


If we fix attention on C, we find that the number а of sequences 
favorable to C is 1, and the probability of termination in C is 
a/y = 1/4. Hence, relative to this category, we have p = 1/4 for 
W-1 апа а = 3/4 for И = 0, with dispersion ра = 3/16. 

Suppose for comparison that we decide to use a weight method for 


this problem as indicated below. 


Q © 
Gece HEHE 
© 


Now sequences of random numbers are limited to length 7 = 1, and 


we have the results 


-198- 


Google 


il 
| 
—. 
no 
DJ 
u 
- 
—. 
no 


О implies А 


ll 
- 
—. 
№ 
Q 
Il 
an 
—. 
no 


1 implies А 


Again fixing attention оп С, 


p, = 1/2 for W = 0 
p, = 1/2 for W, = 1/2 


and 


a=p, М, + p, М, = 1/4 


о 
И 


Py we + Dy We acu = 1/16 < pq = 3/16 


which illustrates in this very trivial example the type of improvement 


discussed in the preceding section. 


(2) Now suppose we have the simple multiplicative system in- 


dicated by the following scheme: 
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Thus a source particle may terminate їп A with probability 1/2 
or, in the contrary case, it may either terminate in B with probability 
1/2 or produce two particles with probability 1/2. The two progeny 
thereupon have equal chances of terminating in A or B. 

Here the maximum length А of random number sequences of O's and 


ц 


l's is 4, у= 2 = 16, and we have the following correspondence between 


random number sequences and terminal weights: 


0000 А=1 Bz0 
0001 А = 1 В=0 
0010 А = 1 В=0 
00) A-1 B=0 
0100 А=1 В= 0 
0101 Az=1 B=0 
0110 А=1 В=0 
01 A-1 B=0 
1000 A=0 Bel 
1001 A=0 В = І 
1010 А = 0 В = 1 
1011 А-0 В=1 
1100 А = 2 В=0 
1101 А = 1 В = 1 
1110 А = 1 Bel 
1111 А = 0 В = 2 


T 


Let us consider only category B. have by inspection the weights 


and corresponding probabilities 
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Ро = 9/16 wW = 0 
p, = 6/16 Му =1 
Po = 1/16 Wo = о, 


The expectation and dispersion аге 


@ 
l 


"im М, = 1/2 


b= р, We ЕЕ 3/8 


Let us compare this procedure with опе which, in the event of 


multiplication, processes one particle of weight 2 rather than two 


particles each of unit weight, as indicated below. 


Google 


We now obtain the results 


000 А =1 В=0 
001 А=1 В=0 
010 А = 1 В= 0 
011 А =1 В=0 
100 А=0 В=1 
101 А=0 В=1 
110 А=2 B=0 
111 А= 0 В= 2 
and for category В, 
ро = 5/8 Wo =0 
р. = 2/8 W-1 
Py = 1/8 М, = 2 


It will be noted that the revision has resulted in increasing the 
dispersion. That this is not an inevitable consequence of employing 
weighted multiplication, the reader may verify by observing that the 
dispersion for category B becomes O in the following treatment of 


the same problem: 


=202= 


Google 
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АРРЕМОТХ 


SUMMARY OF CERTAIN PROBLEMS RUN ON MANIAC I 


Introduction. We conclude this report with a very brief summary 
of some of the Monte Carlo type problems which we have set up and run 
on the MANIAC I at Los Alamos during the past three years. The coding 
for these problems has been done by some extremely rapid and resourceful 
coders at Group T-7 of this laboratory, whom we wish to thank collec- 
tively, and to whom we refer individually in this appendix апа LA-2121. 

It is TE that the reader unfamiliar with the method will find 
in this brief list of problems some indication of the actual situations 
in which the methods of the text find application. It may also perhaps 
help those contemplating use of the method in problems of their own to 
form some judgment of its possibilities. 

The collection included here naturally omits the many classified 


problems which have been completed. 
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Problem 1 


COMPTON COLLISIONS IN A SPHERICAL SHELL 


Coder: В. С. Schrandt 


Requester: G. H. Best and H. C. Hoyt 


15" 


source: Monoenergetic photons directed radially outward from Ro’ 


Medium: Free-electron "gas" in spherical shell of radii Ro < R 


Various initial energies were studied. 
Photon parameters: х, у, Z, ч, у, W, E, в. 
Physical processes: Compton collisions of photons with free electrons. 


Output: Correlated energy and angle distribution of photons escaping Ri 


together with totals lost to core and lost to energy cutoff. Any photon 


impinging on the inner boundary Ro was considered absorbed., 


Remarks: "Forced first collision" device was used. The angular distri- 


bution desired was that of Ch. VIII, Fig. 53. 


Problem 2 


COMPTON COLLISIONS IN A SOLID SPHERE 


Coder: R. G. Schrandt 


Requester: G. H. Best and H. C. Hoyt 


-205- 


Google 


Medium: Free electron gas in solid sphere of radius Ку. 
Source, photon parameters, and physical processes: As іп Problem 1. 
Output: As in Problem 1 with exception of core losses. 

Remarks: It was desired to study the effect of the core on the escape 
distribution by comparison of the results of Problems 1 and 2. Actually, 
by suitable devices in the flow diagram, it was possible to combine the 
two problems into a single code, keeping two sets of terminal category 
registers and following a photon path further for purposes of Problem 2, 


after it had hit the core and been lost for Problem 1. This is a device 


which cuts machine time almost in half. 


Problem 3 


NEUTRONS IN A SPHERICAL SHELL 


Coder: J. М. Kister 

Requester: J. R. Beyster 

Medium: Hollow spherical shell of heavy nuclei. 

Source: Monoenergetic isotropic point source at center of sphere. 
Neutron parameters: В, м, v. 

Physical processes: Inelastic scattering, regarded as terminal, elastic 
scattering according to given lab. system differential cross section, 


attended by no energy loss. 
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Output: Number of neutrons lost to inelastic collision, number of 
neutrons escaping after и collisions, v= O, l, ..., 10, and y 211. 
Remarks: The distribution function к, was tabulated for scattering at 
a cosine a 2 а, for 32 strategically chosen intervals of the cosine 
range. The problem was run as a check on an analytic method of 

H. A. Bethe based on transport cross sections, the results of the 
latter being computed by J. В. Beyster. The Monte Carlo procedure 
was coded in two different ways, one using the weight method for in- 
elastic termination, the other not. The results were in excellent 
agreement and checked Bethe's result. Transmission was obtained 
without forced first collision and agreed closely with the analytic 
value. The problem is written up in LA-1503, A Monte Carlo Deter- 
mination of the Escape Fraction for a Scattering Spherical Shell with 


Central Point Source. 


Problem 4 


ENERGY INDEPENDENT SCATTERING IN A CYLINDER 
Coder: R. L. Bivins 


Requester: M. Walt 


Medium: Solid cylinder of homogeneous material. Problem run for many 


heavy elements. 
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Source: Monoenergetic neutrons in parallel beam incident on lateral 


surface of cylinder in direction normal to axis. 
Neutron parameters: х, у, 2, ч, V, W, И, v. 


Physical processes: Inelastic collision, treated as terminal, elastic 
scattering at source energy governed by differential cross section. No 


energy loss assumed on elastic scattering. 


Output: Transmission, loss to ве collision, number of emergent 
scattered neutrons not hitting detector band, classification of neutrons 
hitting band zone i with v=1, 2, 3, and y 24 collisions. Detector 
band was coaxial with cylinder and calibrated in 5° zones. Geometry 


was that of Fig. 55. 


Remarks: Weights were used for transmission and inelastic losses. 
Elastic scattering probabilities Хх were used as indicated in Problem 3. 
It will be noted that many of the problems included are of this general 
nature, the Monte Carlo results being used in correcting experimentally 
determined differential cross sections and cross sections of other 


processes for multiple collision effects in thick targets. 


Problem 5 


ENERGY DEPENDENT SCATTERING IN A CYLINDER 


Coder: В. L. Bivins 


Requester: М. Walt 
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Medium: Solid cylinder of homogeneous material. Many light elements 


were run. 


Source: Same as Problem 4, 
Neutron parameters: x, y, 2, u, V, W, E, в, W, v. 


Physical processes: Inelastic collision regarded as terminal, elastic 


scattering governed by differential cross sections e (E) in C.M. system. 


Output: Transmission, loss to inelastic collision, loss to energy cutoff, 
emergent neutrons failing to hit detector band (same geometry as in 
Problem 4 and Fig. 55), and those hitting band in zone i with V 2 1, 2, 


and v 2 3 collisions. 


Remarks: Weights used for transmission (forced first collision device) 
and for inelastic losses. Differential elastic scattering handled by 
the von Neumann device using o (eo) = A, +B H+ Chu” ав indicated 
in Ch. V, 85. Tables stored for free path Ag and for coefficients 


Ag? Be? Ce for suitable energy groups g. 


Problem 6 


CYLINDRICAL SHELL WITH PARALLEL SOURCE 


Coder: В. С. Schrandt 
Requester: J. D. Seagrave 


Medium: Cylindrical shell of homogeneous material. Problems run for 


Ср, С, and CH, shells. 


2 2 
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Source: Monoenergetic parallel beam of neutrons incident on lateral 
surface of cylinder in direction normal to axis. Various source energies 


were studied. 


Neutron parameters: x, у, 2, u, V, W, E, в, И, V. 


Physical processes: Elastic scattering on C, D, Н according to C.M. 

differential cross sections e). Such scattering on Н was assumed 
* 

isotropic; C scattering was treated by a fit о (g,p) = A + c at; 


while for D, the double interpolation method of Ch. V, 85 (Fig. 3T) was 


х 
used on the basis of a stored table d^ n° 
3 


Output: Transmission, loss of "degraded" neutrons to energy cutoff, 
emergent undegraded scattered neutrons failing to hit band, and clas- 
sification of scattered undegraded neutrons hitting the band (geometry 


of Fig. 55) into the following categories: 


(a) hits on band zone i with energy above given upper bound є, 
(b) hits on band zone i with energy below given lower bound G 
(с) hits on band zone i with energy between e: and €; with у= 1 


(а) hits on band zone i with energy between ы and €, with и>], 


Problem 7 


14 MEV NEUTRONS IN А CYLINDRICAL SHELL 


Coder: В. С. Schrandt 


Requester: L. Rosen and L. Stewart 
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Medium: Solid cylinder of homogeneous scattering material. Problems run 


for two heavy elements (Bi, Ta). 


Source: Parallel beam of monoenergetic (14 Mev) neutrons incident on 


base of cylinder. Geometry of Fig. 52. 


Neutron parameters: x, у, 2, ц, у, м, и, Е, f (group index for elastic 
scattering differential cross section purposes), g (group index for free 
path and probability of elastic scattering constants), h (group index 


for escape classification). 


Physical processes: At most, two collisions were permitted, escape being 
forced after a second collision, with the direction obtaining as a result 
of that collision. The branching process of Fig. 61 was involved. The 
assumptions made for the types of collisions were as follows: 

(a) Elastic collision: no energy loss, lab. differential cross 


sections e. (a) for five energy ranges f. 


(b) d-process. Energy of emergent neutron equidistributed between 
5 and 12 Mev. Process occurs only at lh Mev. Опе differential 


cross section given for lab. angle &. 


(c) (n-@n) process. Two neutrons emerge with energy distributions 
n, = Е! exp(- к/т and n, = E' ехр(- E'/T,) on OS E! < 5, 
Isotropic emergence in lab. system. Process occurs only at 


14 Mev. (Cf. Ch. V, 815, for method.) 
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(d) Inelastic collision. One neutron emerges, isotropically in 
lab. system, with energy distribution n = E' exp(- E/T) where 
T =k V(E/14). Process occurs only below 12 Mev. (Cf. Ch. V, 


811, for method.) 


\. 
„ v4 e 
«c7 
$75zE"z12 (d 
[/ 
< Sb, 
N IS 
v S (n. 
Ф z 
„$ 
AY 
7 
$ LS 
A Pd 
E - 14 Mev 5 = Е! = 12 (d-process Э 
хр 
о T 
Ф * 
NS 
>) 
Ф 
s Vw 
J к 
p“ 
0 
ts - 
ER 
el, 
Fig. 61 


Output: Energy and angle distribution of escape. 
Remarks: Symmetry of system about source distribution obtains, and 
counter was at infinite distance so that the optimum conditions of 


Fig. 52 obtained. (Cf. discussion of Ch. VIII, 83.) 
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Problem 8 


CYLINDRICAL SHELL WITH POINT SOURCE 


Coder: В. L. Bivins 
Requester: L. Cranberg 
Medium: Cylindrical shell of homogeneous material, Various heavy 


elements run. 


Source: Monoenergetic point source incident on lateral surface of cyl- 
inder from an external point midway between base planes. (Cf. Ch. II, 


854, апа Fig. 11.) 


Neutron parameters: х, у, 2, и, у, м, E (three discrete values Е, Е, 


Ej). 
3) 

Physical processes: А typical element (Bi) had two excited state energy 
levels at А у = .9 and A 


with energy Е colliding with the nucleus can thus scatter elastically, 


р = 1.65 Mev above the normal state. А neutron 


with no energy loss, or inelastically, emerging with energy E - 4, if 


Е = A, or with energy E - 4, if E 2 4,. Since the source energy was 
2.5 Mev, the: histories of Fig. 62 were possible. Inelastic scattering 
was assumed isotropic in the lab. system, while elastic scattering was 


assumed to obey a given differential cross section o(a) applicable at 


all energies involved. 
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Output: At most, two collisions were allowed. Neutrons having a third 


collision were classified in a terminal category L At most, one in- 


3° 
elastic collision was permitted. Neutrons suffering a second inelastic 
collision were classified in a terminal category I5. Thus all escaping 
E,, E 


neutrons emerge with one of three possible energies E Of 


00797 
these, the neutrons hitting the detector band (geometry of Fig. 55) іп 
angular zone i and energy E were tabulated as N, (E) for each Е. 

Also recorded were the transmission, escapes not hitting the band, 


total hitting band after one elastic collision, and the total nitting 


the band after one collision of any kind. 
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Problem 9 


A SCINTILLATION COUNTER 


Coder: R. L. Bivins 
Requester: В. F. Taschek and М. J. Terrell 
Medium: Cylindrical shell containing various compounds of Н, С, О, and Cd. 


Source: Monoenergetic point source at center of cylindrical hole (vacuum) 


with various deterministic directions. 
Neutron parameters: х, у, 2, а, V, м, E, б. 


Physical processes: H-capture, Cd-capture, elastic scattering isotropic 
in C.M. system for Н, О, Cd, and according to a differential cross section 


2 
c =A +С и for C. 
g " " 


Output: Loss from shell surface, H-capture, Cd-capture, and losses to 


energy cutoff classified according to 6 radial and 12 height zones. 


Remarks: The purpose of the problem was the determination of the space 
distribution of degraded neutrons referred to, for use as input in а 
subsequent age theory calculation (not performed by us) connected with 


sensitivity of a scintillation counter design. 
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Problem 10 


A "LONG COUNTER" 


Coder: В. L. Bivins 

Requester: A. W. Schardt 

Medium: Cd cylinder (medium ш = 1, parts р = 1, 2, 4 of Fig. 63) con- 
taining B (gaseous boron compound under pressure) and vacuum space, and 
surrounded by a hydrocarbon (medium m = 2, parts р = 3,5) in cylindrical 


outer container as shown. 


b G OO gle TECHNICAL REPORT ARCHIVE & IMAGE 


source: Monoenergetic neutrons in parallel beam impinging on various 
radial zones (p,p') of base of outer container. Various source energies 


studied. 


Neutron parameters: х, у, Z, и, у, м, В, в, W, m, p. 


Physical processes: H-capture, Cd-capture, B-capture, elastic scattering 


on H, Cd isotropic in C.M. system, C as in Problem 9. 


Output: Escape from Cd-base (p = 1, р = 4), escape from hydrocarbon 
surface (р = 3, 5), H-capture, Cd-capture, B-capture, losses to energy 


and weight cutoff. 


Remarks: The purpose was to estimate the efficiency of B-capture as a 
function of source energy and distance from axis in connection with 
design of a long counter. Weights were used throughout, neutrons 
passing through B cylinder were attenuated exponentially, and the 
multiplication device of Ch. V, §16, was used. This problem probably 
represents the most demanding one statistically that we have set up. 
Enormous samples were involved, and the actual running of the machine 


was done by Schardt and co-workers. 


Problem 11 


A PROBLEM CONNECTED WITH NEUTRINO DETECTION 


Coder: R. L. Bivins 


Requester: Е. Reines and С. Г. Cowan 
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Medium: Plane slab of homogeneous CH, 5 


with Cd-layer at 2 = О, Two different thicknesses of Cd were studied. 


on the z-interval - Н < z € H, 


The slab thickness 2H was essentially infinite. 


Source: Monoenergetic neutrons at various heights 2 and with various 


directions from OZ. А series of source energies was studied. 


Neutron parameters: 2, w, E, T, W, and a parameter indicating the status 


of the neutron relative to two problems run simultaneously. (See Remarks.) 


Physical processes: H-capture at thermal energy (only), Са-сар+иге, 
elastic scattering on H isotropic in C.M. system at energies above 
thermal, and isotropic in lab. system at thermal with no further energy 
loss. Elastic scattering on C was assumed isotropic in the lab. system 


with no energy loss at all energies. 


Output: Classification of Cd-captures in .2 psec intervals from time 
(t = О) of leaving source to time cutoff of 30 usec, 150 intervals in 


all, hydrogen capture, loss to time cutoff. 


Remarks: The following conventions on Cd wall capture were adopted. 

The thicker wall captures all neutrons impinging on it at energies 

<.4 ev, while the thinner wall captures all neutrons impinging at 

energies <.34. Тһе thick and thin Cd-wall problems were run simulta- 
1 t 

12 tt Naso and №, ee Niso 


for the two time distributions and following the path of a neutron 


neously by keeping separate counters N 


further for purposes of the thin wall problem after it had been clas- 


sified terminally as captured in the thick wall problem. 
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Problem 12 


AN @ DETERMINATION 


Coder: J. M. Kister 

Requester: К. В. Lazarus 

Medium: Bare homogeneous Oy sphere. 

Source: Initially ап arbitrary distribution of neutrons in energy 
groups g = 1, 2, 3, radial zones 2 = l, ..., 10, and cosine intervals 


Jj = l, ..., 10, 300 categories in all. 


Neutron parameters: R, w, V - yE, &. 
Physical processes: The usual processes for а fissionable element were 


treated by the transfer matrix method of Ch. V, 812. 


Output: The problem was run in cycles of fixed time length Ат, the 
output Noz, of one cycle serving as the number of neutrons to be 
processed in category (g,z,j) during the next. The problem run was 

close to critical, so that no re-normalization was involved. (cf. Ch. II, 


88, апа Ch. IV, 89.) 
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Problem 13 


NEUTRON FLUX IN AIR 


Coder: В. б. Schrandt 
Requester: J. Hall, R. G. Wagner, and G. M. Wing 
Medium: Sphere of homogeneous air zoned for computational purposes by 


С < 2 
radii Ry Вр < R3 


Source: Monoenergetic isotropic point source at center R = О. Various 


initial energies were studied. 


Neutron parameters: R, w, E, в, W, 2. 


Physical processes: Elastic scattering on O and N assumed isotropic in 


С.М. system. 


Output: Loss to energy cutoff in each zone, loss from system (at R3), 
loss to weight cutoff. These were the only terminal categories. How- 


ever, the purpose of the problem was to determine total flux at R., R 


а? 


В Counters were, therefore, provided to tally all neutrons crossing 


3° 


each boundary ш each energy group g in outward direction and in 


inward direction. (Cf. Ch. IV, 8.) 


Remarks: Weights were used in connection with the device discussed in 
Ch. V, 87, Ғог avoiding loss of trajectories to energy cutoff. Thus no 


trajectory terminated until loss to weight cutoff or loss from Ra. 
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Problem 14 


ESCAPE DISTRIBUTION FROM A CARBON SLAB 


Coders: M. В. Wells and В. б. Schrandt 

Requester: T. B. Taylor 

Medium: Plane slab of С опо & 2 < 2. 

source: Monoenergetic neutrons upwardly directed at 2 = О in cosine 


distribution. 
Neutron parameters: z, w, Е, g". 


Physical processes: Elastic scattering on C governed by a differential 


% 
cross section of form о (g,u) = Ag + c, n. 


Output: 13 energy group classification escapes at z = 0 and 2 = 2. 
Energy cutoff losses classified according to position in five sub- 


intervals. 


Problem 15 


ESCAPE DISTRIBUTION FROM A SPHERICAL SHELL 


Coders: М. B. Wells and В. С. Schrandt 
Requester: С. L. Longmire 


Medium: Spherical shell of hydrocarbon оп В. < К < R 


l 2° 
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Source: Monoenergetic neutrons directed into shell at Ri in cosine 


distribution. 
Neutron parameters: В, w, Е, g". 


Physical processes: Elastic scattering on H isotropic in C.M. system, 


on carbon as in Problem 14. 


Output: 13 energy group classification of escapes at В, escapes at 


Ry» loss to energy cutoff. 


Problem 16 


A ROCKET MOTOR 


Coder: В. L. Bivins 

Requesters: С. I. Bell and С. L. Longmire 

Medium: Cylindrical container of C, U, H, with C of uniform density, 
H and U of prescribed densities in each of a set of shell zones, de- 
fined by 7 radial and 15 height intervals. This cylinder was sur- 
mounted by a coaxial cyclinder of H and Be at specified densities. 
and these two cylinders were surrounded by а cylindrical shell of Be. 
The U-density distribution was subject to change from problem to 


problem, of which & whole series was run. (Cf. Fig. 6h.) 
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Fig. 64 


Source: Isotropic source in fission energy spectrum, spatially distrib- 
uted as prescribed among the 105 zones. The spatial distribution also 
varied with the problem and was so normalized that integral numbers Ng. 


of neutrons were involved in each source. 
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Neutron parameters: х, у, 2, и, м, W, E, 8> ^; W, and & sector parameter 
р= 1, 2, 3. (See Fig. 64.) The zone parameter de really а pair i,j 


which define the associated radial and height zones. 


Physical processes: Be-capture at "thermal" energy; U-capture, fission, 
elastic and inelastic scattering; elastic scattering on C and Be assumed 
isotropic in lab. system, using transport cross sections for free path 
contributions but with correct energy losses. Elastic scattering on H 


isotropic in C.M. system. 


Output: Losses from system classified according to surface of escape; 
Be-capture; U-capture; number ЇЇ! of fissions in each of the 105 spatial 


zones. Fission was regarded as а terminal event. 


Remarks: There was no loss to energy cutoff; neutrons degrading to 
"thermal" energy were allowed to scatter isotropically in the lab. 
system without further loss of energy. А neutron was followed until 
capture, fission, or escape occurred. The purpose of the problem was 
to determine а U-density distribution together with а fission source 
distribution in space which would produce & steady state in the sense 


that the output N' and input N, should satisfy the equations 


7 7 


where P is the fission multiplication constant. Such problems arise 
in the design of nuclear rocket motors, the hydrogen of the problem 


being present in the form of & propellant gas. 
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Problem 17 


PHOTON DIFFUSION IN A REACTOR 


Coder: В. С. Schrandt 
Requesters: М. E. Battat and B. M, Carmichael 
Medium: Cylindrical assembly consisting of inner fuel region, tamper, 


and top and bottom reflectors of given electron densities. 


Source: Prompt У -rays in exponential energy distribution n(E)dE 
= К ехр(- 1.01 E)dE, .2 < E < 5 Mev, isotropic in direction, and with 
given radial and heightwise photon density distributions throughout 


the fuel. 


Photon parameters: x, y, 2, u, м, м, E, g, m (medium index indicating 


2 ue 


region of occupancy: fuel, tamper, or reflectors), R' , К", H', Н" 


defining radial and height boundaries of region m. 
Physical processes: Compton scattering of photons on electrons with 
energy losses due to such scattering and to cutoffs (varying with m) 


chosen at the point where the photoelectric effect becomes an 


effective absorption. 


Output: Distribution in energy апа angle of photons escaping each of 
outer surfaces, energy deposition within the system in 82 spatial 


regions. 
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Problem 18 


А THICK Zr TARGET PROBLEM 


Coder: В. С. Schrandt 

Requesters: L. P. Stewart and L. P. Rosen 

Medium: Solid Zr cylinder. 

Source: Parallel beam of monoenergetic (14 Mev) neutrons incident on 


base of cylinder. 
Neutron parameters: х, у, 2, ч, у, м, v, E, f, g, h (cf. Problem T). 


Physical processes: Those of Problem 7 except that the (-2n) process 
was replaced by & more general process involving production cross 
sections for neutrons in each of two energy distributions with a 


total expectation of 2.32 neutrons. 
Output: As in Problem 7. 


Remarks: The problem was handled by a revision of the code of the 


problem cited. 


Problem 19 


A THICK C TARGET PROBLEM 


Coder: В. б. Schrandt 


Requesters: L. P. Stewart and L. P. Rosen 
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Medium: Solid carbon cylinder. 


Source: Monoenergetic (14.1 Mev) neutrons in parallel beam incident on 


base of cylinder. 


Neutron parameters: u, v, м, x, y, Z, E, в, v. 


Physical processes: Elastic scattering according to given differential 
cross sections curves on four ranges between 14 and O Mev; capture; 
inelastic scattering for levels 4.4, 7.6, 9.6, and 11.2 Mev (the latter 
an average for levels >9.6), inelastic scattering at 4.4 level for 
incident energy >10 Mev governed by given c (E,0), all other inelastic 


scattering assumed isotropic in C.M. system. 


Output: Distribution in 28 energy groups апа 10? angular intervals 
(with vertical) of neutrons escaping after one or two collisions, 


total number of second collisions, 


Remarks: First collision forced, neutrons forced out after a second 
collision. Method of Ch. V, 89, used to determine parameters after 


inelastic collision. 


Problem 20 


HEAVY WATER EXPERIMENT 


Coder: В. С. Schrandt 


Requester: J. М. Grundl 
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Medium: Heavy water in (А) spherical shell, (В) hemispherical shell. 
Both problems run simultaneously. A series of runs involved various 


inner and outer radii. 


source: Isotropic point source at center of system in fission-neutron 


energy distribution, .l © ES 10 Mev. 


Neutron parameters: u, у, W, х, у, Z, Е, g, M (parameter indicating 


trajectory lost relative to problem (B) or not). 


Physical processes: Elastic scattering on D and О assumed isotropic in 
C.M. system; transport cross section based on this assumption used for 


free path. 


Output: Рог problems (A) and (B), loss to outer surfaces, loss to energy 
cutoff, classification in 25 energy groups of neutrons re-entering central 


cavity with normal angles in each of four angular ranges. 
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